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Abstract: We report on a calculation of the cross-section for Higgs boson production in 
gluon fusion in association with a hadronic jet at next-to-next-to-leading order (NNLO) in 
perturbative QCD. The computational technique is discussed in detail. We show explicitly 
how to employ known soft and collinear limits of scattering amplitudes to construct subtrac- 
tion terms for NNLO computations. Cancellation of singularities is demonstrated numerically 
for the collinearly-subtracted gg ^ H + j cross-section through NNLO and the finite dgg^Hj 
cross-section is computed through 0{a^) as a function of the center-of-mass collision energy. 
We present numerical results for the gluon-fusion contribution to Higgs production in asso- 
ciation with a jet at the LHC. The NNLO QCD corrections significantly reduce the residual 
scale dependence of the cross-section. The computational method that we describe in this 
paper is applicable to the calculation of NNLO QCD corrections to any other 2—7-2 process 
at a hadron collider without modification. 
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1 Introduction 

The ATLAS and CMS experiments at the Large Hadron Collider (LHC) have discovered a 
new particle with a mass of approximately 125 GeV [1, 2] whose properties are consistent 
with that of the Standard Model Higgs boson. Continuing studies at the LHC are focusing on 
the detailed understanding of the quantum numbers of this particle [3] and its couplings to 
gauge bosons and fermions (see e.g. [4] and references therein). The successful completion of 
this task is crucial for determining if the new particle is indeed the long-awaited Higgs boson 
or instead some other state. 

A reliable understanding of coupling constants cannot occur without accurate theoretical 
predictions for the main Higgs boson production and decay processes. Arriving at such 
predictions requires the computation of higher-order QCD corrections, since they are known 
to affect Higgs production rates and decay branching fractions in a significant way. In fact, 
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for gluon-initiated processes gg ^ H + X , where X is a state with zero, one or two hard 
jets, the next-to-leading order (NLO) QCD radiative corrections are known to be so large [5- 
12] that next-to-next-to-leading order (NNLO) QCD computations are important for reliable 
phenomenology. 

Unfortunately, current computational technology only allows NNLO QCD computations 
for the case of Higgs boson production in association with zero jets [13-17]. Extending this 
result to one or more jets will lead to a refined analysis of the pp ^ H ^ W~^W~ process, 
since in that case final states with different jet multiplicities are treated as different processes 
in order to optimize search strategies. The information about the relative significance of the 
Higgs boson production in association with zero, one or two jets is currently extracted from 
available fixed-order computations supplemented with resummations of the most important 
terms in the perturbative expansion [18-21]. Explicit NNLO computations for multi-jet 
processes will be indispensable for understanding the reliability of these predictions. 

The other motivation for this work is of a more theoretical nature. One can argue that 
the framework of perturbative QCD that has been developed since the late 1970's has proven 
to be one of the most important areas of particle physics phenomenology. Indeed, it is im- 
possible to imagine contemporary high-energy physics without hadron collider physics whose 
proper description is intimately related with parton shower event generators, sophisticated 
fits of parton distribution functions, fixed-order perturbative calculations and the like. Much 
of our understanding of perturbative QCD is based on how soft and collinear singularities 
cancel in suitable "infra-red safe" quantities, since this defines short-distance observables that 
can be calculated in perturbation theory. Currently, there exists an interesting gap in this 
understanding. On one hand, general theorems [22, 23] ensure that this cancellation occurs 
in suitably defined quantities to all orders in perturbation theory. On the other hand, we only 
know how to use those ideas for generic computations of infra-red-safe observables at leading 
and next-to-leading order in perturbative QCD [24, 25]. It is still not entirely clear how to 
construct a general calculational scheme for two- and higher-loop computations. 

It is important to point out that, in spite of the fact that a generic computational scheme 
is not available, a fairly large number of NNLO computations for various processes have 
already been performed [16, 17, 26-39], but until very recently such computations always 
utilized a particular aspect of a specific process. Such aspects included a small number of 
final-state particles, or their color neutrality, or absence of color-charged particles in the initial 
state, or even the fact that all matter particles in a particular process were massive. A generic 
algorithm that is valid irrespective of the details of the process under consideration was not 
worked out. This situation is somewhat peculiar, because following the successful development 
of generic methods for NLO computations [24, 25], it was generally felt that the development 
of similar methods for NNLO computations would be relatively straightforward. For this 
reason, about ten years ago many authors calculated infra-red and collinear limits of generic 
QCD amplitudes [40-47] that are potentially relevant for NNLO computations, and a large 
number of two-loop 2 — )• 2 scattering amplitudes became available [48, 49]. Unfortunately, 
since it proved harder than expected to develop a working scheme for NNLO computations. 
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these infra-red and collinear limits were never used for their intended purpose. 

Recently, important steps towards developing a general computational scheme valid 
through NNLO QCD that, at least in principle, is applicable to processes of arbitrary mul- 
tiplicity, were made by Czakon [55, 56], who suggested to combine the ideas of sector de- 
composition [57-59] and Frixione-Kunszt-Signer (FKS) [25] phase-space partitioning. These 
results were used to obtain NNLO QCD corrections to the cross-sections for qq — )• tt [60, 61] 
and qg — ?■ tt [62] processes. A similar computational scheme was also developed in Ref. [63] 
where it was applied to the calculation of NNLO QED corrections to Z — )• e~^e~. We note 
that parallel developments in the antennae subtraction technique [64] have recently led to the 
calculation of the NNLO QCD corrections to the leading-color all-gluon contribution to di-jet 
production at the LHC [65]. 

While the above results mark important progress in developing a suitable framework 
for NNLO computations, a large number of technical aspects still needs to be understood 
and worked out. It is best to do this by considering a realistic example with all the ensuing 
complications. This is the purpose of the current paper. We consider the hadro-production of 
the Higgs boson in association with one hadronic jet at NNLO in QCD. To make calculations 
as simple as possible, we work in pure gluodynamics, i.e. QCD without light fermions. We 
note that an understanding of how to compute the NNLO QCD corrections to H+j production 
is instructive because this process possesses all non-trivial aspects of a generic NNLO QCD 
problem. Indeed, 

• there are colored particles in the initial state; 

• there are colored particles in the final state; 

• already at leading order, the total cross-section for this process does not exist unless a 
jet algorithm is specified; 

• this process exhibits the most general structure of infra-red and collinear singularities, 
since these singularities occur due to radiation of gluons in the initial and final states; 

• singular collinear splittings g ^ gg involve non-trivial spin correlations; 

• the number of Feynman diagrams that we need to compute is large. ^ 

The only "non-generic" feature of the process that we consider is the high symmetry of 
the final state which simplifies the bookkeeping and speeds up the computation. We feel, 
however, that having this simplification is useful in the first step in the development of the 
new technology and that it does not affect the generality of the method that we describe in 
this paper. 

The remainder of the paper is organized as follows. In the next Section we describe the 
setup of the calculation. In Section 3 we discuss the parametrizations of the phase-space 

^For example, the process gg Hggg is described at leading order by 230 diagrams while the gg — > Hgg 
process at one-loop is described by 603 diagrams. 
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for leading, next-to-leading and next-to-next-to-leading order computations. In Section 4 we 
explain how singular limits of amplitudes are used. In Section 5 we describe how 0{e) parts 
of the relevant amplitudes can be computed using helicity methods. In Section 6 we describe 
the numerical implementation of our method. In Section 7 we discuss some tests and show 
the results of the computation. We conclude in Section 8. Some useful formulae are given in 
the Appendix. 



2 The setup 

We are interested in the computation of NNLO QCD corrections to the process g + g ^ 
H -\- g, where the Higgs boson can decay into arbitrary particles. To compute this and 
related processes, we use the QCD Lagrangian, supplemented with a dimension-five non- 
renormalizable operator that describes the interaction of the Higgs boson with gluons in the 
limit of very large top quark mass 

^ = -^Gj'jG^'^)'^'^ - A^^,,ifG(fjG('^)''^'^ (2.1) 

Here, G^j^u is the field-strength tensor of the gluon field and H is the Higgs boson field. 

Matrix elements computed with the Lagrangian of Eq. (2.1) need to be renormalized; to 
do so, two renormalization constants are required. The first one relates bare and renormalized 
QCD coupling constants 

«<"V?5. ^ ^■a.Z.^. ^ I - f (^) . (f - I) il^y . .... ,2.2, 

Here, Og is the strong coupling constant in the MS scheme evaluated at the renormalization 
scale Se = {^TT)~'^e~'^^ , 7 = 0.5772 is the Euler constant and 

,. = 11^, A = ,2.3, 

o 

are one- and two-loop contributions to the QCD /3-function computed with the Lagrangian 
of Eq. (2.1). A'^c = 3 is the number of colors. We note that Eq. (2.3) is only valid in a theory 
without light fermions, as defined by the Lagrangian Eq. (2.1). 

The second renormalization constant ensures that matrix elements of the HGG dimension- 
five operator are finite. It reads 

In the above formula, G{as) is the Wilson coefficient of the HGG operator in the MS scheme 
[51] 



^ , 11 /as\ /as\2 
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+ 0{as,Nj), (2.5) 
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where mt is the mass of the top quark. We emphasize again that the displayed result is only 
valid in the approximation when no light fermions are present in the theory. 

Renormalization of the strong coupling constant and of the effective Higgs-gluon coupling 
removes ultraviolet divergences from the matrix elements. The remaining divergences are of 
infra-red origin. To remove them, we must both define and compute infra-red safe observables, 
and absorb remaining collinear singularities by renormalizing parton distribution functions. 
We now discuss these two issues. 

Generic infra-red safe observables are defined using jet algorithms. For the calculation 
described in this paper we employ the /c^-algorithm. This algorithm belongs to the class 
of sequential jet algorithms. It requires specification of the minimal transverse momentum 
of the reconstructed jets p±^j and the minimal "angular" distance between two partons, 
^Rij = \/ {Vi - VjY + i'^i - V'i)^) where y = l/21n(S + pz)/{E - p^) is the rapidity and 
ip is the azimuthal angle of a parton. Once AR and p±j are specified, the jet algorithm 
maps a set of parton momenta onto a set of jet momenta in such a way that jet momenta 
are stable against soft and collinear parton splittings. The Kinoshita-Lee-Naunberg theorem 
[22, 23] then ensures that observables constructed from jet four- momenta are determined by 
short-distance physics and can therefore be computed in QCD perturbation theory. However, 
because massless colored partons are present in the initial state of the partonic process gg — >• 
H + X, the infra-red and collinear cancellation is not complete, even in the presence of a jet 
algorithm. Collinear singularities associated with gluon radiation by incoming partons must 
be removed by additional renormalization of parton distribution functions. We describe how 
to perform this renormalization in what follows. For definiteness, we focus our discussion on 
the production cross-section of a Higgs boson and a jet in pure gluodynamics. 

We denote the UV-renormalized partonic cross-section for the production of the Higgs 
boson and a jet in a gluon fusion by a{xi,X2), and the collinear-renormalized partonic cross- 
section by a{xi,X2). Once we know a{xi,X2), we can compute the hadronic cross-sections by 
integrating the product of a and gluon distribution functions over xi and X2 




(2.6) 



The relation between a and a is given by the following formula ^ 



a = T ^ (g) (g) r 




where the convolution sign stands for 



1 









^We show this relation for pure gluodynamics; if more species of partons are present, Eq. (2.7) becomes a 
matrix equation. 
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The collinear counter-terms are defined as 



Os\ „ (a 



r = *a--)-gjr. + g) r„ (2,9) 

The relevant sphtting functions and their convolutions are given in the Appendix. We write 
the UV-renormalized partonic cross-section through NNLO as 



2 



with 



and the collinear-renormalized partonic cross-section as 

We note that the collinear-renormalized cross-section is finite. We then use Eq. (2.7) to obtain 



+ r(i + .) + r(i + e)' 

(2) _ -(2) _ ra ® _ ® Ta _ Ti ® a(o) ® Ti Ti (g) ^(i) a(^) ^ Ti 

r(i + e)2 r(i + e)2 r(i + e)2 + r(i + e) + r(i + e) ■ 

Although finite, the a^*^ still depend on unphysical renormalization and factorization 
scales because of the truncation of the perturbative expansion. In the following, we will 
consider for simplicity the case of equal renormalization and factorization scales , ^Ur = ^/ = /i. 
The residual [i dependence is easily determined by solving the renormalization group equation 
order- by-order in a^- The equation reads 



The //-derivative of the right hand side can be computed using the known evolution equations 
for the strong coupling constant and the gluon density 



Solving these renormalization group equations, we get 

+ Pi? « Pi? « + ofl P(? P(? + 2P(? ^ ^ Pi?) , 



(2.15) 



(2.16) 
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where = a^^\as{fJ-), fJ.) and L12 = In/^f/^l- 

It follows from Eqs. (2.13,2.16) that, in order to obtain cj^^^ at a generic scale, apart 
from lower-order results we need to know the NNLO renormalized cross-section a^'^^ and 



convolutions of NLO and LO cross-sections with various splitting functions. Up to terms 
induced by the renormalization, there are three contributions to a^'^^ that are required: 

• the two-loop virtual corrections to gg — )• Hg; 

• the one- loop virtual corrections to gg ^ H + gg; 

• the double-real contribution gg ^ H + ggg. 

We note that helicity amplitudes for all of these processes are available in the literature. 
The two-loop amplitudes for gg — )■ Hg were recently computed in Ref. [50]. The one-loop 
corrections to gg — )■ Hgg [52] and the tree amplitudes for gg — )■ Hggg [53] are known. 
Moreover, in the two latter cases, these amplitudes are available in the form of a Fortran code 
in the program MCFM [66]. In principle, they can be just taken from MCFM and used with 
no modification in another numerical program. 

Since the above discussion implies that all ingredients for the NNLO computation of 
gg ^ H + jet are available and, in fact, have been available for some time, it is important to 
understand what has prevented the community from performing this and similar calculations. 
In fact, the main difficulties with NNLO calculations appear when we attempt to combine the 
different contributions, since integration over phase-space introduces additional singularities 
if the required number of jets is lower than the parton multiplicity. To perform the phase- 
space integration, we must first isolate singularities in tree- and loop amplitudes. It required 
a long time to establish a convenient way to do this. 

The computational method that we will explain shortly is based on the idea that relevant 
singularities can be isolated using appropriate parametrizations of phase-space and expansions 
in plus-distributions [25, 59]. To illustrate this point, we consider the integral 



1 




(2.17) 







where the function Fix) has a well-defined limit lim Fix) = Fii)). We would like to construct 

x->-0 

the Laurent expansion of / in e. This can be accomplished by writing 




(2.18) 



so that 



1 




F(0) F{x) - F{0) 



ae X 



— ae 



F{x) - F{0) 



X 



ln(x) + ... 



) 



(2.19) 
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The above equation provides the required Laurent expansion of the integral /(e). We note 
that each term in such an expansion can be calculated independently from other terms. 

To use this approach for computing NNLO QCD corrections, we need to map the relevant 
phase-space to a unit hypercube in such a way that extraction of singularities is straightfor- 
ward. It is intuitively clear that correct variables to use are the re-scaled energies of unre- 
solved partons and the relative angles between two unresolved (collinear) partons. However, 
the problem is that different partons become unresolved in different parts of the phase-space. 
It is not immediately clear how to switch between different sets of coordinates and cover the 
full phase-space. 

We note that for NLO QCD computations, this problem was solved in Ref. [25], where 
it was explained that the full phase-space can be partitioned into sectors in such a way that 
in each sector only one parton (i) can produce a soft singularity and only one pair of partons 
(ij) can produce a collinear singularity. In each sector, the proper variables are the energy 
of the parton i and the relative angle between partons i and j. Once the partitioning of 
the phase-space is established and proper variables are chosen for each sector, we can use an 
expansion in plus-distributions to construct relevant subtraction terms for each sector. With 
the subtraction terms in place, the Laurent expansion of cross-sections in e can be constructed, 
and each term in such an expansion can be integrated over the phase-space independently. 
Therefore, partitioning of the phase-space into suitable sectors and proper parametrization 
of the phase-space in each of these sectors are the two crucial elements needed to extend this 
method to NNLO. In the next Section we discuss these issues in detail. 



3 Phase-space parametrizations and sector decomposition 
3.1 Phase-space for leading order processes 

We now discuss how to parametrize the leading-order phase-space for the process gi + g2 ^ 
H + g-^. This will be needed both for the leading-order cross section and for the NLO virtual 
and NNLO double virtual corrections, so it must be computed in d = 4 — 2e dimensions. We 
note that the integration over the leading order phase-space is not singular, because of the 
requirement that a jet is observed. We work in the center-of-mass frame of the two incoming 
gluons, so that their momenta are parametrized as 

pi = ^(1,0,0,1), P2 = ^ (1,0, 0,-1). (3.1) 

The center-of-mass collision energy is denoted by ^/s and the mass of the Higgs boson is 
denoted by mn- The production cross-section, averaged over spins and colors of the two 
colliding gluons, is written as 

dogg^H+g = -^^dUY>Si2-,m\Mgg^gH\'^ X Fj, (3.2) 

where Fj is the "measurement function" that restricts the integration to the region of phase- 
space where there is an identified jet. The amplitude Aigg^gH describes production of an 
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on-shell Higgs boson in hadronic collisions and dLipS]^2^H3 is the Lorentz-invariant phase- 
space. 

The parametrization of the phase-space dLipsi2-s.3_H' in Eq. (3.2) is easily obtained by 
integrating over the momentum of the on-shell Higgs boson and then over the center-of-mass 
energy of the gluon §3. We find 



dO^'^"^VlH dcos 63 



dLipSi2_^3^ = — 8(2^)^-2 ~[}~~^)' P±,H = Em^^smOs, (3.3) 

where i?max = {s — mjj) / {2^/s) . With this parametrization, the momentum of the gluon (73 
reads 

P3 = -Emax (1,^3) , (3.4) 

where 723 = (sin ^3 cos (^3, sin ^3 sin (/?3, cos ^3), and the Higgs boson momentum is obtained 
from momentum conservation: pn = Pi + P2 — Ps- Note that in Eq. (3.3), p_l^h is the 
transverse momentum of the Higgs boson relative to the collision axis. 

Before proceeding further, we note that the azimuthal angle ips of the emitted gluon is 
a dummy variable, since neither matrix element squared nor the measurement function Fj 
depend on it for our choice of pi and p2- Hence, we can rotate it away by taking the gs 
momentum to be 

P3 = -E'max (1, sin 6*3,0, cos 6*3) , (3.5) 

and integrate over dQ^^ in Eq. (3.3). Once this is done, we can set this solid angle to its 
four-dimensional expression to simplify calculations at higher orders. This is legitimate to do 
as long as we can identify this angle, associated with global rotations of final states in the 
plane transverse to the collision axis, when parametrizing higher-multiplicity phase-spaces. 
Nevertheless, to maintain sufficiently general leading order kinematics, we re-introduce the 
azimuthal angle ip3 and keep it to generate momenta of the gluon and the Higgs boson. 
This amounts to writing dQ^ — )■ dips in Eq. (3.3) and then using Eq. (3.4) for the gluon 
g3 momentum. In addition, since the transverse momentum of the Higgs boson is an ob- 
servable quantity, the diff'erential cross-section da/dp±^H should be finite for each value of 
P±,H to all orders in perturbation theory. Hence, we can divide the cross-section by p^% 
without changing the final result. This amounts to removing this factor from the phase- 
space parametrization at both leading and higher orders in perturbation theory. Putting all 
these remarks together, we conclude that we can choose the leading order phase-space to be 
"four-dimensional," 

.T- dcos6l3d(/?3 / / mjj\ dxsdx 4 , 

dLips,2_3H ^ 32^2 \'-—) = \^-—) (^-^^ 

where we introduced cos ^3 = 1 — 2x3 and 933 = 27rx4 to parametrize the momentum of gluon 
g3 as given by Eq. (3.4). We must remember to normalize NLO and NNLO phase-spaces to 
P^H consistency. 
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Although we will not discuss this in any detail in this paper, we note that it is straight- 
forward to include decays of the Higgs boson. Indeed, because the Higgs boson momentum 
is an observable quantity, all singularities should cancel out in the differential cross-section 
da/dpH- Once this differential cross-section is known and because the Higgs boson is a scalar 
particle, so that no spin correlations are present, we can easily turn da/dpH into quantities 
such as dcr/dp!yjdp!y2 by letting the Higgs boson decay in its rest frame and then boosting the 
four-momenta of the two photons into the center-of-mass frame. 

3.2 Phase-space for next-to-leading order processes 

In this Section, we consider the parametrization of the phase-space for the process gi + (72 — ^ 
H + g3 + 54. This process represents a real-emission contribution to the production cross- 
section of the Higgs and one jet at next-to-leading order. It is also important for the NNLO 
computation where integration of one-loop corrections to gg — )■ Hgg amplitudes over the ggH 
phase-space is required. 

As we already explained in the Introduction, a good parametrization of the phase-space 
dLipSg^g,^_^jjg,^g^ shouM f aclht at c the extraction of singularities from the matrix elements of 
the process 5152 — ^ Hg^g^. Of particular importance are collinear singularities. We can 
compare two cases: i) 54 is emitted collinear to gi] ii) 54 is emitted collinear to 53. In the 
first case, it is easiest to extract the singularity if the z-axis is chosen to coincide with the 
direction of the gluon gi and in the second case with the direction of the gluon (73. This 
immediately tells us that a suitable parametrization of the phase-space should depend on the 
kinematics of the process. As we mentioned in the Introduction, this is the main idea behind 
the FKS subtraction method [25]. 

Following Ref. [25], we note that the first step towards a convenient phase-space para- 
metrization is the phase-space partitioning. The goal of such a partitioning is to create sectors 
where one and only one gluon or one and only one pair of gluons can become unresolved. 
Once we know which gluon or which pair of gluons can produce singularities, we choose 
the energy of the potentially soft gluon and the relative angle between the two potentially 
collinear gluons as the primary variables for the phase-space parametrization in the given 
sector. To illustrate this procedure, we begin by removing the symmetry between the two 
gluons in the final state by separating them into "resolved" and "unresolved" ones. To this 
end, we introduce the following function of transverse momenta of gluons 5^3 and (74, 



= ^^'^ , J/^, (3.7) 



and write 



— dLipsi2-^34H = 7TTdLipsi2_>34H (a('^^)+A(32) ^dLips Ag. (3.8) 

In the last step we used the fact that the phase space, the kinematic constraints on final-state 
particles and all matrix elements are symmetric with respect to permutations of gluons (73 and 
(74. Given the structure of the damping factor Ap^, it is clear that singularities of the matrix 
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element related to gluon are unimportant, and we only need to consider cases when gluon 
g'4 becomes either soft or collinear to one of the three hard directions defined by the momenta 
gi,g2 and (73. Note that g^ and 54 cannot both be soft, or collinear to the collision axis at 
the same time, because we require a jet in the final state. To separate the collinear-singular 
regions, we introduce another partition of unity 

1 = Af ) + Af ) + \ (3.9) 

In Eq. (3.9), we use 

~ 14 24 I 14 S4 I 24 S4 ' J) ^ 7^ ^) 4, (3.10) 

where p*-' = 1 — nj • ftj and fti is the three-vector that parametrizes momentum direction of 
the particle i. Again, the A^'^*^ are labeled in such a way that the subscript indicates a pair 
of particles that can become collinear without forcing the angular damping factor to vanish. 
Inserting this partition of unity Eq. (3.9) into the phase-space of Eq. (3.8), we obtain 

-dLipSi2^34// ^ XI ^^^P^12i34if 

dLipsit'ig^^ = dLipsi2^34^A(^>f \ (3.11) 

i=l 

The above decomposition defines pre-sectors that we will refer to as Sc'-'^'^ A phase-space 
parametrization for each of these pre-sectors is chosen in such a way that the soft and collinear 
singularities that are relevant for that pre-sector can be extracted in the easiest possible way. 

We now describe these parametrizations explicitly. In general, we will parametrize the 
phase-spaces by splitting them into "regular" and "singular" parts 

dLipsg'^g^^ = A(4>f )dLipsQ(i2)^3H X [d54](^*). (3.12) 

The regular NLO phase-space is the same for all pre-sectors. It includes all particles except 
the (potentially soft) gluon g^. We write it as 

where we have introduced the notation Q = pi + P2 — Vi and = -Egg (1, ^3). Also, 

-Egg = — ^ " , ns = (sin 6*3 cos 933 , sin 6*3 sin (^3 , cos 6*3 ) , 

2(Qo-Q-"'3) (3.14) 

COS03 = 1 — 2x4, sin 03 = — cos^ ^3, (^3 = 2-71x5. 

Following the discussion of the leading order phase-space parametrization, we have dropped 
the e-dependent part of the integral over azimuthal angle of the gluon g^, and have normalized 
the remaining e-dependent part of the phase-space to the transverse momentum of the Higgs 
boson. 
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Parametrization of the singular phase-space depends on the pre-sector. To explain this, 
we begin by considering pre-sector Sc'-^^^ To parametrize the singular phase-space for this 
pre-sector, we note that, thanks to the damping factors, the singularities occur when 54 is 
collinear to gi or when §4^ is soft. Hence, it is convenient to choose the parametrization where 
the energy of the gluon and the relative angle between the three-momenta of 54 and gi are 
basic variables. The azimuthal angle of the gluon 54 is conveniently defined relative to the 
plane formed by the gi and 53 three-momenta. We therefore write 

Pi = Eg^ (1, sin 6*4 cos 994, sin 64 sin (^4, cos 64) , (3.15) 

where (/?4 = (f4 + cp^. The singular phase-space reads 



E^~MEg,d COS g4dy4(sin2 ^4 sin^ y4)-^dl7('^-3) 
) = . (3.16) 



The {d — 3)-dimensional solid angle does not enter any of the scalar products and therefore 
can be integrated away. We write 

Eg4^ = Em!,xXi, COS 6*4 = 1 - 2x2, (3.17) 

where .Emax is introduced after Eq. (3.3). The singular phase-space for Sc*^^^^ becomes 

1(41) _ 2-2^i^f~^^ i_2e-ef-. ^ . d(^4(sin^(v?4 " V^s))"' . .^ 

[d54j^ ' - -femax (-27r)'^-^ ^ ^ (1-^2) — • (3.18) 

/ d(^4(sin2((^4))-^ 



We use 

2 ,„ \-e _ r,l-2e - ' 1 1 





and write ip^ = 2ttx3 to find 

^2 ^4 



d<^4(sin^ ^4)-' = 2^-''B ( - - e, - - e ) , (3.19) 



3 

X x}"2'x^'(l - X2)"'(sin^((/?4 - V^s))"' n 



(3.20) 



1 

Combining everything, we find the expression for the phase-space of the pre-sector Sc^^^^ to 
be 

dLipsiS34H = Norm x PS^ x pg-.^^i^^^^Jad^ ^ ^^2^^] ^ (3 2i) 



where PS = IGE^^^Sf sin^ 6*3(1 - X2){sm'^{ip4 - 953))/^!,^ and 

r(l + e) / 7r2 2 3 ttS 

Norm = / Z>n 1 e - 2C3e + — e 

(47r)^/2 V 3 ^'^ 90 

pc _ -^max-^3 ^(41). (4) 

87rH^s-E4{l-n,.n4))'^' 



(3.22) 
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These equations allow us to generate four- momenta of all final-state particles. Indeed, 
a set of random numbers Xi,..,x^ gives us momenta of the gluon and the direction of 
the unit vector that parametrizes the momentum direction of a "hard" gluon (73 in the 
center-of-mass frame of colliding gluons. Using this information, we can find the energy of 
the gluon 53 and determine the momentum of the Higgs boson from momentum conservation. 

The phase-space of the second pre-sector Sc^^^^ dLips^g^g^^^^ is parametrized in a similar 
way, except that we now need a simple parametrization of the relative angle between gluons (72 
and 54. Therefore, we write cos 84 = —1 + 2x4. This is the only change that occurs at the level 
of momentum generation and everything else, including the phase-space parametrization, can 
be borrowed from Eq. (3.21). 

The phase-space parametrization for the third pre-sector Sc^^'^^ requires some changes. 
The main difference with respect to the previous cases is that now the collinear direction 
corresponds to the "hard" final state gluon (73, which means that we need to choose the 
relative angle between and g^ as the primary variable for phase-space parametrization. In 
the reference frame where the momentum of gluon g^ is along the z-axis, the direction of the 
gluon (74 is chosen to be 

^4,3||z = (sin04 cos (^4, sin ^4 sin(/?4, cos ^4) . (3.23) 

The phase-space parametrization employs angles 64 and ip^. The momentum of the gluon 54 
in the center-of-mass reference frame is obtained by rotating Eq. (3.23) in the x — z plane 
by ^3 and in the x — y plane by (f^. We parametrize the energy of the gluon (74 and its 
relative angle with respect to (73 using Eq. (3.17). We conclude that the parametrization 
of Sc^^'^-' phase-space coincides with Eq. (3.21) except that in PS~^, we should substitute 
(sin2(y?4 - (^3))-^ ^ (sin2(994))-^ 

The above formulae can be used to construct phase-space parametrizations for next- 
to-leading computations or for the calculation of the one-loop corrections to gg ^ H + gg 
process. In the latter case, one should be careful since it is customary for one-loop virtual 
corrections to be normalized with the factor 

r(i + 6)r(i-6)^ 

= (4vr)2-r(l - 2e) ' ^''''^ 
If we choose the normalization in such a way that one power 

of r(l-Fe)/(47r)'^/2 ig factored out 
per loop, the expression for Norm in Eq. (3.22) changes. To use Eq. (3.21) for the computation 
of real- virtual corrections, we should make the following replacement there 

Norm ^ NormRv = cpNorm = ^ (l - - AC^e^ + ^e^) . (3.25) 



3.3 Phase-space for next-to-next-to-leading order processes 

In this Section we consider the partonic process 91+52 H + g^ + g^ + g^ and discuss 
how to generate the phase-space in a way that facilitates the extraction of singularities. We 
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begin with a discussion of the phase-space partitioning. Similar to the one-loop case, we first 
partition the phase-space in a way that allows us to identify the "hard" gluon by writing 

4f = ' ^^i^^' i,j,k€[3A,5]. (3.26) 

Because A^'^^^ + A^^^) + A^^^^ = 1, we can use this partition of unity and the symmetry of the 
phase-space, the measurement functions and the matrix elements with respect to permutations 
of gluons (73 , 54 and 55 , to write 

idLipsi2_^345 = |dLipsi2^^345 (a^^) + A^^) + A^f) = IdLipSi^^^g^^ A^^^^) 

In the last step we introduced the energy ordering of the two gluons; this allows us to remove 
the final symmetry factor. 

We must next partition the phase-space to extract collinear singularities. To do so, we 
closely follow the discussion of the next-to-leading order case in the previous Section. We split 
the phase-space into nine different sectors that we denote by the possible collinear directions 
of the gluons 4 and 5. We have three triple-collinear sectors 4||5||i, with i = 1,2,3 and six 
double-collinear sectors 4||z(g)5||j, where j & [1,2,3]. To write the weight for each of the 
nine sectors, we introduce the auxiliary quantities 

3 3 

di€[4,5] = Pij ^ die[4,5]k = ^ Pij , d^^ij = + p^i + p^j . (3.28) 

Denoting the weight of a sector where gluon 4 is allowed to become collinear to gluon i and 
gluon 5 to gluon j by wa-^j, we write (/c / n / 4 7^ 5 7^ f 7^ j) 



I _ PikPinP^kPbn 
W4i^3\i=j- ^^^^ 



dik din J V^5fc d^n 



- + —) (— + —)^+(- + —) (± + ±)^ (3.29) 

dii d^k J \d5k d^n J d^^ni \dAi din J \d^k d^n J di^ki 



^dik din J \d^i d^kj di^in \dik d^nj yd^i d^nj di5ik 

and (A; / n / 4 / 5 / i, / / m / 4 / 5 / i) 

I P4kP4nP5lP57n f , \ f , ^ \ P45 Qn\ 

Using Eq. (3.27), we decompose the phase-space as 

^dLipSi2_,/f345 = X]'^^^P^12^i^345' (3-31) 
ae5 
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where S = [(41; 51), (42; 52), (43; 53), (41; 52), (42; 51), (41; 53), (43; 51), (42; 53), (43; 52)] and 

dLipsS2^^345 = dLipSi2^//345 ^px^ ^(^94 " ^95) ^a- (3.32) 

We now discuss the parametrization of the phase-spaces for individual pre-sectors. Be- 
cause of the Ap^^ factor, we consider gluon as part of the regular phase-space and gluons 
94) 95 3-5 part of the singular phase-space. Regular phase-spaces are the same for all pre-sectors 
and are parametrized in the same way as at NLO in Eq. (3.13), except that the vector Q in 
that equation becomes Q = pi + P2 — Pi — P5- 

We begin with the triple-collinear sectors. We have three such sectors Sc^^*'^*\ i G [1, 2, 3]. 
In these sectors, singularities can appear if gluons (74^5 are soft, and if they are collinear to the 
direction 1%, or to each other. The phase-space parametrization should enable us to extract 
all of these singularities. We will start the discussion with the triple-collinear initial sector 

g^(41;51)_ 

The first step is to find independent degrees of freedom, which is non-trivial because we 
have to perform computations in dimensional regularization. To illustrate this point, we use 
d-dimensional rotational invariance to choose the momenta of five gluons as follows 

Pi,2 = ^(l,0,0,±l;0), 

P3 = Eg^{l, sin 6*3 cos if 3, , sin 6*3 sin ips , cos 6*3 ; 0) , (-3 33) 

Pi = Eg^ (1, sin 64, 0, cos 64; 0) , 

P5 = Eg^ ( 1 , sin ^5 cos , sin 6^ sin tp^ cos a , cos 65 ; sin 6^ sin tp^ sin a) . 

Note that these momenta are shown as five- dimensional vectors; the fifth component corre- 
sponds to one of the axes in the {d — 4)-dimensional space. The angle a parametrizes leakage 
into the {d — 4)-dimensional vector space. Note also that we have chosen to give the {d — 4)- 
dimensional component to the softer of the two gluons. The reason for this choice will be 
explained shortly. With this parametrization, the angular part of the phase-space becomes 

dJ)g-i)dJ^g-^)dJ7(f ^) ~ d[cos^3](sin2 03)-^d^3(sin2(^3)-^dO(f ^) 

X d[cos^4](sin2 04)~'d!^J'f"2) d[cos ^sKsin^ ^5)"' (3-34) 
X d(^5(sin2 if5)-'d[cosa]{sm^ ar^-'dnf-'^l 

We can generalize the momentum parametrization in Eq. (3.33) by rotating all momenta 
in the xy-plane by the angle 9J4. Obviously, the momenta of the incoming gluons pi^2 do not 
change, while the other momenta become 

Pi = (1, sin 64 cos tpi, sin 6*4 cos (/?4, cos ^4; 0) , 

P5 = E's (1, sin 05 coSq(v34 -|- 995), sin 65 siua^^i + ^5), cos 9^; sin 6^ sin ip^ sina) , (3.35) 
P3 = £'3 ( 1 , sin 63 cos (fs sin 6^ sin 993 , cos ^3 ; 0) . 
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In Eq. (3.35), we have introduced the notation 

coSq,((/94 + ip^) = COS ipi COS — sin 994 simp^ cos a, 

(3.36) 

sinQ,((^4 + (^5) = sin 994 cos 9J5 + cos 994 sin 995 cos a, 

and (/93 = (^3 + (/94 . Note that the phase-space is written in terms of (^3 , the relative azimuthal 
angle of (74 and g^, and that 



Before we express the phase-space parametrization in terms of suitable variables, we make 
a few general comments. We note that our choice of the phase-space parametrization and 
assignment of extra-dimensional components is restricted by two requirements: 

• extra-dimensional components and angles should not complicate the extraction of sin- 
gular limits; 

• extra-dimensional momenta components should not appear in the non-singular matrix 
elements and kinematic constraints. 

It turns out that the parametrization of the momenta in Eq. (3.35) satisfies the first require- 
ment for the triple-collinear sector Sc'-^^'^^^ This happens because the parametrization is 
chosen in such a way that the scalar products pi ■ P4,pi ■ p^^Pi • P5 that can potentially lead 
to singularities in this sector do not depend on the extra-dimensional angle a. 

We now discuss how to satisfy the second requirement. We note that full parametrization 
of Eq. (3.35) is not needed for the highest multiplicity gg — )■ Hggg hard matrix element. 
Indeed, a configuration where all the three final-state gluons are resolved is non-singular, 
hence we can use a d = 4 phase-space parametrization to describe it. We will see explicitly 
below that this amounts to setting a = in Eq. (3.35). Therefore, we only have to explain how 
to satisfy the second requirement in configurations where one or both of 34, g^ are unresolved. 
To this end, we note that in all soft limits this requirement is automatically satisfied. Indeed, 
since Eg^ — ^ implies -E^g — t- 0, in any of the soft limits the gluon momentum with the e- 
dimensional component is not present in the hard matrix element and in kinematic constraints. 
The a-dependence will therefore reside solely in the unresolved phase-space and in eikonal 
factors and splitting functions. It is important that this dependence on a is non-singular, so 
that the numerical integration can be performed in a straightforward way. 

The collinear limits are more complicated. If p^ is collinear to either pi or p^, then 
(p5 = or 65 = 0, which implies that the e-dimensional components of momenta and the 
dependence on a disappear from the matrix elements. On the other hand, this does not 
mean that collinear limits are independent of a. Indeed, such a dependence is present in the 
spin-correlation part of the splitting functions. We must account for that in the computation. 
This can be done in a straightforward way since this dependence is non-singular. Finally, 
consider the kinematic situation where p4 is collinear to pi and p^ is resolved . In this case, 
the matrix element squared becomes 




(3.37) 



\M\' 



1 



{2CAgs)P^n{PA,I^A)M^'{pM,P2,P3,P5)M*'''{pM,P2,P3,P5) 



(3.38) 



Pi - Pi 
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where pi4 = Pi — Pa and K4 is the spin-correlation vector that tehs us how the colhnear 
direction is approached (see Section 4 or [24] for details). Eq. (3.38) implies that the matrix 
element depends on the four- vector p^ and, according to Eq. (3.35), ^5 has e-dimensional 
components. This dependence is unfortunate, since it becomes unclear how to use four- 
dimensional methods, such as spinor-helicity techniques, to simplify calculations of scattering 
amplitudes in that situation. However, when pi||p4 we are left with only three different 
directions ni, ns, n^. We can use d— dimensional rotational invariance to remove any e- 
dimensional components from the matrix elements in Eq. (3.38). To do so, we first remove the 
y-component of P3 by rotating all momenta in the xy-plane by the angle — fs- This rotation 
does not change ~ (1,0,0,1) and p2 ~ (1,0,0,-1). We then perform another rotation 
in the ye-plane, to remove the e-dependent component of the vector p^. Because none of the 
momenta in the matrix element has both y- and e-dimensional components, such a rotation 
does not change pi4,P2 and p^, while it makes p^ four-dimensional. We note that, although we 
rotated away the e-dimensional components of the resolved four-vectors that are used in the 
hard matrix elements, these vectors still depend on the e-dimensional angle a. In addition, 
because of spin correlations, we also must rotate the vector = (0, cos 934, sin 994, 0, 0) that 
enters Pgg in Eq. (3.38). This rotated vector receives e-dimensional components and becomes 
a-dependent. The purpose of the rotation therefore is to move the e-dimensional components 
from the resolved momenta in the matrix element to the splitting function, where it is easy 
to account for them explicitly. Finally, we stress that the very possibility to rotate away 
the e-dimensional components of particle momenta is connected to the rotational invariance 
of spin-summed scattering amplitudes squared in d-dimensional space-time. This seems to 
suggest that the easiest framework in which to implement this techniques is conventional 
dimensional regularization, where the momenta of all external particles and their polarization 
vectors are treated as d-dimensional. We will discuss this point in more detail shortly. 

We now discuss the explicit parametrizations of the relevant phase-spaces. For the sector 
gc(4i;5i)^ the singular phase-space reads 

[d54] [dg5]e{Eg^ - Eg^) = 24+2e(27r)2d-2 '^'^^ [sin^(v'4 - V^s)] ' d cos a [sin^ a] 

X [66]^-" [mil - m)r [mil - r,,)r [A(i - A)]-/2-e (^•^9) 

X (2E„ax)'"'' ^(6 - 6)^ (Uax - 6) d6d6d??4dr?5dA. 

The variables introduced in the above formula parametrize the energies and angles of the 
(potentially) unresolved gluons in the following way 



-^54,95 — -E'maxCl,2, Cmax — min 



1, 



l-il 



1 - (1 - m2/s)^i7745 



(3.40) 



and 



\m-m? .2 .^\m-m? 

r/45 = , sm (f^ = 4A(1 - A) 



D = riA^m- 2??4% + 2(2A - i)\Jmmii - m){i - %)■ 
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The two variables 774^5 are scalar products of the reference direction vector ni and the vectors 
that parametrize directions of the two gluons 



2r/4,5 = 1 - "^4,5 • ni. (3.42) 

The parametrization of triple-collinear phase-spaces in Eq. (3.39) is still too complicated to 
extract all singularities; further decomposition is required. This is achieved by a sequence of 
variable changes that we describe below, following Refs. [55, 56]. Specifically, we split the 
triple-collinear initial-initial sector into five sectors 

dLips(4^^5^)=^dLips(^^^5i'^). (3.43) 

i 

To project onto individual contributions, we need to perform the following changes of variables 



Sc(4i;5i,i) . 


6 


= Xi, 


6 


= XiXYas.xX2j ^4 = 


X3, V5 


X3X4 
2 ' 




Sc(4i;5i,2) . 




= Xl, 




= XiXmax2;2) ^4 = 


X3, V5 


= ^3 (1 - 


X4\ 

2 / 


Sc(4i;5i,3) . 




= Xl, 


6 


= XiXt^is,xX2X4:, 7?4 


_ X3X4 

2 


, V5 = X3 




Sc(4i;5i,4) . 




= Xl, 


6 


= XiXt^s_xX2i ^4 = 


X3X4X2 
2 


, m = X3 




Sc(4i;5i,5) . 


6 


= Xl, 


^2 


= XiXt^s_xX2i ^4 = 


^3(1- 


X4\ 

yj, % 


= 2:3 



(3.44) 



We also write A = sin^(7rx5/2). This change of variables introduces a factor of vr in the 
normalization of the phase-space that is included in the expressions below. The integration 
region for X5 is always between zero and one. 

We also note that the {d — 4)-dimensional angle a introduces singularities in the phase- 
space parametrization. To take care of them, we calculate the integral over this angle, 

1 1 

/dcosa 1 r dxg r(— e)^ 

[sin2a]-i+^ "21^7 xl+'{l-xgy+^ " 2i+2^r(-2e) ' ^ ' 







and write 



d[cosa] _ r(l — 2e) dxc 



sm a 



2a]i+^ " 2r(l-e)2^ x^+^(l-X9)i+^ 



r(l-2e) ^ dx9(l-X9)- 



(3.46) 



where cosq = 1 — 2x9 and in the last step we used the symmetry of the matrix element 
with respect to X9 -H- 1 — xg, to simplify the integrand. We can expand Eq. (3.46) in plus- 
distributions. Such an expansion does not introduce additional poles in e. We find 



Xg 



5(X9) 



1 

Xg 



+ .. (3.47) 

+ 
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Note that the first term in the expansion corresponds to a = 0, which reduces the parametriza- 
tion of momenta of all final-state particles to their four-dimensional limits. The "extra- 
dimensional" momenta components and the "extra-dimensional" angles appear with an ad- 
ditional suppression in e, but because of infra-red singularities, they start contributing to 
differential cross-sections already at 0{e~'^). 

For each of the five sectors Sc^^^'^^'*\ we write the phase-space in the form 



/_ X 9 ^ dx 
dLipsg.51 ~ Norm X PS^,,PS-^ x ^ H ^^'^ ^ 11 ^ 



9 fc=5 j=lx. ' 



b^^ bi''> 6<') ftW' 

-| O/O li/Q 



(3.48) 



(i) (i) 

Below we present the functions PS^^j, PSj and the exponents 0^=1 4 and bjj^-^ ^ for each of 
the sectors. First, we note that the normalization factor is common to all sectors; it reads 



Norm 



' ' - - 2C{3)e' + ) • (3-49) 



r(i + 6) 

(4vr)d/2 

We also note that we can write 
PS„ = J:,^^PS„, PS. = 12?M^^!^#M(i^.in^,,,3)P5., (3.50) 

where (p^s = (p4, — (fs- The expressions for the exponents and the phase-space factors for each 
of the five sectors read (we suppress the sector label everywhere in the equations below) 

Sector Sc^^^'^^'^) : {ai = 4, 02 = 2, 03 = 2, 04 = 1}, {61 = 4, 62 = 2, 63 = 2, 64 = 1}; 
PS (1 - f )^Lx 



2iVi(x3,f ,A)' 



„2 



_ ^max (1 - ^) A(l -X){1-^)'{1- X3) 
2iV2(x3,f,A) 

Sector Sc(^^'^^'2) : {01 = 4,02 = 2,03 = 2,04 = 2}, {61 = 4,62 = 2,63 = 2,64 = 2}; 
^2 

pg "^max 



4iVi(rE3,l-f ,A) 



PS 



X 



2 

max 



(l-X3)(l-f) (l-X3(l-f))A(l-A) 



4iv2(x3,l-f ,A) 



Sector Sc^^^'^^'^^ : {01=4,02 = 2,03 = 2,04 = 8}, {61=4,62 = 2,63 = 2,64 = 8}; 



2 (-1 _ X4\ 

'^^'"-2iVi(x3,f,A)' 

_ xLx(l - ^3) (1 - ^) (1 - f )' A(l - A) 
2iV2(x3,f,A) 
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Sgg^ Sp41;51,4 . {0^=4^^2=3,03 = 2,04 = 1}, {61=4,62 = 3,63=2,64 = 1}; 

2 (-1 _ X2X4 \ 
pg _ -^maxV-^ 2 



pe _ ^max(l - ^) (l - (1 - X3) A(l - A) 



gp(41;51,5) . {0^=4^^2 = 2,03 = 2,04 = 2}, {61=4,62 = 2,63 = 2,64 = 2}; 
_ ^2 

PS. 



4iVi(a;3,l-f ,A)' 



_ _ xLx(l - ^3)(1 - f )(1 - X3(l - f ))A(1 - A) 



4iv2(x3,l-f,A) 

(3.51) 



The function A^^i reads 



iVi(x3,X4,A) = 1 + X4(l - 2x3) - 2(1 - 2A) Vx4(l - X3)(l - X3X4). (3.52) 

The above phase-space parametrization is such that the Hmits Xj — )■ 0, i = 1...4 of the matrix 
element squared can be easily computed; we will discuss this in more detail in the next 
Section. In the remainder of this Section, we will focus on the phase-space parametrization 
of the other pre-sectors. 

We note that the phase-space parametrization for the triple-collinear pre-sector Sc^^'^^ 
is constructed in exact analogy to Sc^^'^^. The only difference is that the collinear direction 
is now n2 = (0, 0, —1) instead of ni = (0, 0, 1). This means that, in terms of the 77- variables, 
angles of gluons (74^5 relative to the collision axis are given by cos ^4^5 = — 1 + 2^74^5. 

The construction of the phase-space parametrization for the triple-collinear pre-sector 
gp(43;53) slightly more involved, since the collinear direction now is the direction of the 
gluon (73. It is therefore convenient to write momenta of 54 and 55 in the reference frame 
where 53 is along the z-axis. We write 

pi''> = Eg, {1,0,0,1;0) 

^4^^ = £'34 (1, sin ^4 COS 994, sin ^4 sin 994, cos ^4; 0) , (3.53) 
p^^^ = Eg^ ( 1, sin 6*5 coSq, ((/?45 ), sin eissinQ, (9945), cos 6*5 ; sine's sinews sin a) , 

where 9945 = 994-1-995. In this sector, the scalar products whose vanishing leads to singularities 
are Ps • Pi,P3 ■ P5 and Pi • P5. It is easy to see from Eq. (3.53) that these scalar products are 
independent of a. The phase-space for Sc^^'^^ depends on two relative angles 994 and 995, 
so that Lips43.53 ~ (sin^ 994 sin^ 995) To get the momenta in the center-of-mass frame, 
we rotate these vectors first in the xz plane by ^3, and then in the xy plane by 993. The 



- 20 - 



parametrization of the singular phase-space is similar to what we have discussed in connection 
with Sc^^^'^^^, except that the collinear direction now is n^. 

We finally turn to the discussion of the double-collinear sectors. First, consider the sectors 
where collinear singularities arise from emission along two incoming particles, (Sc^'^^'^^) and 
gg(42,5i)-j^ In such sectors, scalar products whose vanishing may create singularities are p4,5 -pi 
and P4.5 ■p2- Vanishing of the scalar product -p^ cannot lead to singularities in this sector, 
see Eq. (3.30). With this in mind, we parametrize momenta of the three final-state gluons in 
the center-of-mass frame as 

P3 = Eg.^ (1, sin 6*3 cos (^3, sin 6'3 sin (^3, cos 6*3; 0) , 

Pi = Eg^ (1, sin 6*4 cos((^3 + (^4), sin 6^ sin((^3 + (^4), cos 64; 0) , (3.54) 
P5 = Eg^ (l,sin6'5COSa(v33 + (^5),sin6'4sina(¥53 + (^5), cos 6^5; sin 6*5 sin 995 sin a) . 

The phase-space is parametrized in terms of the relative angles (^4 and (^5. We find 
[dg,][dgMEg, - Eg.J = ^^-^^f-'^^f-'^e{i, - ^2)0{U.. - 6) 



X d6d6^J~^'^2"^'d cos Bid cos 65 (sin^ e^yism^ 65) 
d(/?4(sin^ ((^4))"'^ d(^5(sin^ (f^)"'^ d [cos a] 



(3.55) 



j dv94(sm^v'4) ^ J d(^5(sm^(^5) ^ 


where Eg^^g^ = -EmaxCi,2- We now change variables £,1 = xi,^2 = xiX2Xma,x, cos 64^5 = 1-2x3,4, 
(^4,5 = 27rx5,6 ™d cos a = 1 — 2xq. We use symmetry with respect to xg — )• 1 — xg to simplify 
the expression for the phase-space. We obtain 

dLips4i.52 ~ Norm x PS^PS"' x JJ ^^'^ ^ 11 "tSt ^ [^i'3;2'x33x4*j . (3.56) 

^9 k=5 j=l ^ 



The normalization factors read'^ 



Norm 



r(i + 6) 

(4vr)d/2 



^ / 7r2 37r4 
1 - -e2 - 2C(3)63 + —e' 
2 ^ 40 



PS. = ^ /^^|- PS., (3.57) 

and the exponents read 

{01=4,02 = 2,03 = 1,04 = 1}, {61 = 4,62 = 2,63 = 1,^4 = 1}. (3.58) 



^We note that in Ref. [56] the double coUinear sectors are further spht by an additional partitioning of 
energy and angle variables. We find that such a partitioning is unnecessary. 
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The other type of double-colhnear sectors that need to be considered is the initial-final 
one. We focus for definiteness on Sc^^^'^^^. The momenta read 

P3 = Eg^{l, sin 6^3 cos tps , sin 6^3 sin 993 , cos 6*3 ; 0) , 

P4 = Eg^ (1, sin ^4 cos ((^3 + (^4), sin 6*4 sin ((^3 + (^4).cos04) , (3.59) 

Pg^'* = Eg^ (1, sin ^5 cos (^5, sin 6^ sin (^5 cos a, cos 6*5; sin 6^ sin (^5 sin a) . 

Note that and p4 are given in the center-of-mass frame, while is written in the reference 
frame where is along the z-axis. To obtain p^ in the center-of-mass frame, we rotate it by 
63 in the xz-plane and by tps in the xy-p\ane. The phase-space is identical to Eq. (3.57). The 
discussion of all other double collinear sectors proceeds along the same lines. 

4 Singular limits 

In this Section, we describe the extraction of singular limits. We begin with the next-to- 
leading order computation. We note that we will not discuss the most general case from the 
point of view of color correlations; instead we will make use of the fact that we are studying 
Higgs boson production in association with a jet and so the number of colored particles never 
exceeds five. This feature leads to simplification of the color correlations in soft limits. We 
will make use of these simplifications in what follows. 

4.1 Limits at next-to-leading order 

Consider, for definiteness, the NLO sector Sc'-'^'^^ . The phase-space for this sector, dLips^2^34jy, 
is given by an expression similar to Eq. (3.21), where X2 parametrizes the relative angle be- 
tween 54 and 53. We have to integrate the matrix element squared \A4gg^Hgg\'^ over the 
phase-space. The integration has the form 

1 

l+L l+e ■■■ ^ F{xi,X2, ...), F{xi,X2,...) =[xlx2]\Mgg^Hgg?', (4.1) 
^ ^ 

where the ellipses denote the measurement function, regular parts of the phase-space, various 
damping factors and possible additional arguments of the function F. All of these things 
are not important for discussing the structure of singularities which is shown explicitly in 
Eq. (4.1). The singularities correspond to xi = or X2 = 0, and the function F[xi,X2, ■■) 
is finite in those limits. The integral in Eq. (4.1) is calculated using an expansion in plus- 
distributions, as we explained in Section 2. It follows that in order to perform the integration 
in Eq. (4.1), we need to understand values of the function F(xi, X2, ..) in cases when one (or 
both) of the two first arguments vanishes. 

Consider first the xi = limit. According to the phase-space parametrization described 
in Section 3.2, xi = implies that 54 is soft: Eg^ = 0. In the soft limit, the matrix element 
is written as a product of a reduced matrix element and the eikonal factor 

\Mg,g,^Hg,gf « 2C (l^i,^ + I'S;^ + ^^S^ \Mg,g,^Hg,\\ (4.2) 
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^i = ^..(P^) = 7-^\7^> (4.3) 



where 

{Pi-Pk){Pj -Pk) 

is the eikonal factor. To calculate the soft xi — ?■ limit, we note that the eikonal factor is 
quadratic in p4 = (1,714) ~ xi and so it is easy to compute the required limit. We obtain 



F(0,X2, ...) = — — + — + — \Mgg-^Hg^\ , (4.4) 
Anax VP14P24 Pl4 P24/ 

where we traded X2 for P34/2 which is valid in sector Sc^^'^-'. We note that potential singular- 
ities that correspond to gluon 54 being collinear to gluons gi or §2 are apparent in Eq. (4.4); 
these singularities are removed by the angular damping factor Eq. (3.10) for this sector. 

The second singular limit we have to consider is X2 = 0. In sector Sc^^^^, X2 = means 
that gluon 54 is collinear to gluon (73. The corresponding limit reads 

\Mgg^Hgg\' - f^^if ^(-,^)-^,",^H,-^;r-.H,, (4.5) 

where z = Eg^/{Eg._^ + Eg^) and 

P%a\z,Ki,e) = -g^, + + 2(1 - e)z(l - z)k^^^k^^, (4.6) 

is the gluon splitting function. The vector ^4^^ is the normalized remnant of the momentum p4 
that parametrizes the projection of on the plane transverse to the collinear direction which 
in this case is fixed to be the momentum of gluon g^. Because of the chosen parametrization 
of P4 at next-to-leading order, has only four-dimensional components. 

We will now show how to simplify Eq. (4.5). The idea is to trade the sum over the Lorentz 
indices /i and u for a sum over helicity indices. This is achieved by inserting the completeness 
relation ^ 

where g^'^ denotes the metric tensor of the d-dimensional vector space and h is an auxiliary 
vector such that ■ h ^ 0. Next, we write 



^t,u {Z, ^)J'^^gg^Hg-^^^gg^Hg " I «^A«^A . ~ I 



(4.8) 



where the last step follows from the transversality of the physical amplitude M^pz,fi = and 
from K4 • Pa = 0. Repeating the same procedure with the index we find 



P%3\z,e)M^^M*^^ = Y,Pxx'\^,^WxMl„ (4.9) 
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where the sum over physical hehcities in d-dimensional space-time is performed. 

We now explain how to compute P^y^ ■ First, we note that the polarization vectors 
of a gluon with four-dimensional momenta embedded in a d-dimensional space-time can be 
chosen in the following way. We take the polarization vectors to be either four-dimensional 
vectors that describe states of plus and minus helicity, or (d — 4) dimensional vectors of the 
type = (0, 0, 0, 0; 0, .., 1, 0, .., 0), where projection on a single extra-dimensional direction is 
non- vanishing. The helicity-dependent splitting function P^y\z,e) reads 

Pil^\z, e) = Jav [jZ-, + ^) + - - ^)(eA • ^4) {e*y ■ k,). (4.10) 

For regular it polarizations, both terms in Eq. (4.10) are in general non- vanishing. For extra- 
dimensional polarizations, • K4 = 0, because K4 in this case is a four-dimensional vector, 
and P^xi^ ~ for these polarizations. 

We can now compute P^^^\z, e) and calculate the X2 = limit of the function F2{xi,X2)- 
The final result reads 

F2ix,A.:) = ^f^E^iv^(^'^)-^.'«H-^rtH, (4.11) 

where §3 means that the matrix element should be computed with the momentum of the 
final state gluon given by the sum of the momenta of the gluons (73 and 54. We note that 
Eq. (4.11) requires the computation of scattering amplitudes for the gg — )• Hg process when A 
parametrizes an extra-dimensional polarization vector. We explain how to do this in Section 5. 

Finally, we discuss how the vector K4 is computed. This vector parametrizes how the 
coUinear limit is approached in the plane transverse to the collinear direction. For this 
reason, it depends on the considered sector. To make this explicit, we consider the sector 
Sc*-^^^ and write p4 = xp^ + yp^, + k±K4, where ps = (£"3, —ps), K4 • pa = and ^4-^3 = 
0. A simple computation gives K4 = (cos 63 cos ips cos tp^ — sin (^3 sin (^4, cos 6*3 sin (^3 cos 994 + 
cos (^3 sin (^54, — sin 03 cos (^94). The analogous vectors for the other sectors are much simpler. 
For example, for sectors Sc^^^^ and Sc^^'^\ we find K4 = (cos 9J4 , sin (/74 , 0) . We note that these 
vectors are uniquely determined for each of the phase-space points; this allows us to construct 
the correct splitting function and perform the local subtraction of singularities. The quality 
of the subtraction terms so constructed will be studied in Section 7. 

These are the only two limits that are required for a NLO computation. An expression 
for i^(0, 0, ...) can be easily obtained from the soft limit Eq. (4.4), which is non-singular for 
r/34 — >• 0. Note that the collinear limit has a well-known 1/(1 — z) singularity as gluon g^ 
becomes soft; therefore, to compute -F(0,0, ...) from the collinear limit, one has to cancel xi 
in the numerator in Eq.(4.11) with 1/(1 — z) ~ 1/xi in the splitting function. 

4.2 Limits of double-real emission processes 

In this Section, we briefly discuss the singular limits of the double-real emission processes. As 
already pointed out, the phase-space partitioning splits the phase-space into double-singular 
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and triple-singular sectors. Collinear singularities of the double-collinear sectors are given by 
products of gluon splitting functions, because the two unresolved gluons must be collinear 
to different directions. On the contrary, in the triple-collinear sectors, the 1—7-3 gluon 
splitting functions [54] are required to describe collinear limits. For both double-singular and 
triple-singular sectors, soft singularities originate from both double-soft and single-soft limits. 

We begin by discussing the double-soft limit of the gig2 — )• Hg^g^g^ scattering amplitude. 
It occurs when the momenta of 54 and 55 become vanishingly small. In general, double-soft 
limits involve color-correlated matrix elements, but in our case this does not occur. The 
reason is that, once gluons 54 and 55 decouple, the matrix element depends on three colored 
particles, gi,g2,g3. If, following Ref. [24], we denote the color charge of a particle i by the 
operator Tj, color conservation implies 

4+4+4=0. (4.12) 



In addition, the squares of the color charge operators are equal to the Casimir operators of 

9i 



the SU{3) gauge group. For gluons, this means T^. = Ca- Using these two equations, we find 



■ Tg^ = Tg^ ■ Tg^ = Tg^ ' Tg.^ = -Ca/2, (4.13) 

so that all color correlations are absent. As a result, we can use a simple formula for the 
double soft limit {Sp = [12, 13, 23]) 



yij&Sp I \knGSp 



3 

ijeSp i=l 



(4.14) 



We note that Sij{pk) is given in Eq. (4.3) and Sij{p/^,p^) can be found in Ref. [45]. Using the 
parametrization of the NNLO phase-space and the explicit dependence of the momenta on the 
singular variable xi that controls the double-soft limit, it is straightforward to show that all 
the singularities can be resolved in triple-collinear sectors. In turn, this implies that in such 
sectors we can compute any limit of the form F(0, 0:2, X3, X4, ...) from the double soft-limit, 
with no need to further distinguish the Xj=2,.,4 = case from the Xi=2,..,4 / one. 

Another new element at NNLO is the triple-collinear limit. Similar to double-collinear 
limits, the triple-collinear limits are described by the corresponding splitting functions. For 
example, in the case of the final-state triple-collinear splitting when the momenta of all final- 
state gluons become parallel, we find 

\Mg^g2^Hg3g4.g^\'^ » -^M^^^''^'^-^9i92-^^^9345 "^9192^^3345 ' ^"^'^^^ 
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where S345 = (^3+^4+155)^, and (7345 denotes a gluon with momentum ^345 =^3+^4+^5. 
tio: 

(1-6) 



The sphtting function P^^^^"*^^^ was computed in Ref. [54]; it reads 



(519293) 



^•'12 
S123 



4 S12 



Z3 




2(l-Z3) + 4zJ 1-2^3(1-^3 



1 



Z3 



+ 



Sl23(l 



S12S13 



1-2^3 



~|~ /lo /Cq 



Zl{l - Zl] 
1-2Z2 



■S123 ( 4Z2^3 + 2zi(l - Zl) - I 

+ 2(1-6)^ V (l-^2)(l-^3) 



~|~ ( key 



.,L.\ ^2^2(1 



In Eq. (4.16), 



V 23(1 
{Pi+Pj) 



2'^2Z3(1-Z3)" ^ '^2(1-^2 

_ l-2zi(l-zi) 

2=22:3 

+ 5 permutations 



(4.16) 



Z2 



Z3) 



{Pi +Pj +Pkf, and 



t 



ij,k 



Zi -\- Zj 



+ 



Zi -\- Zj 



The relevant vectors in this formula are computed in the following way. For final-state triple- 
collinear splitting, the momentum of the resolved gluon ^3 defines the collinear direction. The 
energy fractions Zi are obtained from energy ratios Zi = Eg^ / Eg where Es = Eg^ + Eg^ + E'gg . 
Similar to the NLO case, for each phase-space point, we compute directions in the plane 
transverse to the collinear direction p-^ along which collinear limits for 94,(75 are taken. We 
denote such directions as ^4,5, respectively. The vectors that enter the triple-collinear splitting 
function read ^ 



E-^K, = ^4(1 



\ u a 
Zi)K2 - Z4Z5K'^, 

Z5)Kg - Z5Z4K2, 



(4.17) 



E~'^k^ 



-Z3Z4k'^ 



z^z^k!^. 



1. 



5 . 5 

It is easy to check that ^ kg. =0, thanks to the energy-conservation condition ^ Zi 

i=3 1=3 

Finally, we note that for the initial triple-collinear limits that correspond to gluons and g^ 
being collinear to incoming gluons gi or g2, the above formulas are valid up to a replacement 



^93 



-Eg^ or Eg^ 



-E. 



92- 



We also require the triple-collinear splitting function in the strongly ordered configuration, 
Sij <C Sijk <^ 1. In principle, we can obtain it by directly taking the limit of Eq. (4.16). 
However, it is also easy to compute it directly. Indeed, in this case the full triple-collinear 
splitting function factorizes into a (spin-correlated) product of ordinary splitting 



P, 



(919293) 



We give these vectors for physical labels of the three final state gluons. 
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functions. We find that in the strongly-ordered S35 ^ S345 ^ 1 limit, the full matrix element 
can be written as 

\M |2 ~ 0(339495) A^*.!^ (4.^R^ 

The strong-ordered splitting function reads 



S.O., (93,94,35) A 1 y 



+ (1 - e)2z5(l - 2;5)(k4 • K5) 

1 — Z4 



1 — Z4 Z4 / \ 1 — Z5 2:5 

K^K4(1 - e)2z4(l - 2:4) X 



, — + ^ + ^5(1 - ^5)) + ^'^^(l - 6)2Z5(1 - .5)^1 



(4.19) 



where Z5 = Eg^/{Eg^ + Eg^), Z4 is defined as before Z4 = Eg^/Eg and ^4^5 are the spin- 
correlation vectors Hi = ZiPs + k_i_Ki + UiPs- 

4.3 Real-virtual corrections 

In this Section we consider the computation of one-loop corrections to the real-emission 
process gig2 — ?• Hg^g^. We will refer to this contribution as "real-virtual". To calculate 
this contribution, we must integrate the interference of the one-loop and the tree-level matrix 
elements for gig2 — s- H + gsg^ over the NLO phase-space. The NLO phase-space was discussed 
in Section 3.2, where we showed how to partition it in such a way that soft and collinear 
singularities can be extracted. Following Section 3.2, we denote the resolved final-state gluon 
as ^3 and the potentially unresolved final state gluon as g4. For each sector, we denote 
the product of the phase-space parameters and the interference of tree- and one-loop matrix 
elements as 

Frv(xi,X2, ...) = x?X2 2Re (M^^^^M^'^^f^^^^) , (4.20) 



3l32-!>-ff9394 9192->J^9394 ^ 

where xi parametrizes the energy of g^ and X2 = (1 — cos^)/2 parametrizes the cosine of 
the angle between the direction of the gluon (^4 and the collinear direction. This direction 
is sector-dependent, and is given explicitly later. The ellipses in Eq. (4.20) stand for other 
parameters that are needed to fully describe the final-state kinematics. 

We must integrate the function Frv over the phase-space of the softer gluon; schemati- 
cally, the integral takes the form 

1 

dxidx2 ~ , -. 

i+2e i+, -^Rv(a^i,X2,..). (4.21) 
^2 

We note that the extraction of singular limits would have been no different from the NLO case 
discussed in Section 4.1, if not for the fact that the function Fry is not well-defined for xi = 
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and X2 = 0. This happens because F^{xi,X2, ■■) contains branch cuts in the hmits xi — )• 
and X2 — )• 0. To make us of the expansion in plus-distributions, we must isolate these branch 
cuts before extracting the singularities. We can accomplish this by writing Fy^{xi,X2) as the 
sum of three terms 



Frv(2;i,X2, 



Fi{xuX2,...) + {xix2) ' F2{xi,X2,...)+X^'^'F^{xi,X2,...), (4.22) 



where the functions Fi{xi,X2, • • • ) are free from branch-cut singularities so that their values 
at xi = or 2:2 = can be computed. To justify the decomposition of Eq. (4.22), we consider 
the limit when the energy of the gluon (74 becomes small. In this limit, the matrix element 
squared for gig2 ^ H + g^g^ factorizes as [67] 



1-^9192-^-^9394 1^ ^ 5s^^^2C^ (2^12,4 +2^13,4 +2^23,4) l^gig-z^^Hgsl^ , 



where the soft factors 2^ij,4 read 



^^JA=TS\ + 2g!^?^CAC^Il;:, + .... 



r(l) 



The function X^-j ^ is given in Eq. (4.3). The one- loop function X^-^-^ reads 

1 r2(l -e)r2(l + e) 



^(1) 



(4.23) 



(4.24) 



(4.25) 



e2r(l - 2e)r(l + 2e) 

where the eikonal factor Sij can be found in Eq. (4.3). 

We next expand Eq. (4.23) through first order in the strong coupling constant to obtain 



2Re(A^«(5i,..,94)A^(°^'*(5i,-,94 
4:gsCrCA 



94-^0 



^(0) .^(0) ^(0) 
-^12,4 "'"-'-13,4 "'"-^23,4 



2Re LM 



9i92-s>^^93 91 



M 



(0),* / ^^ 

9i92^-f/93V5l) ••) ff3j J 



(4.26) 



+ 



7-(i) I j-y^j I 7-' 



'12,4 



(1) 



13,4 



^23,4 



Using the explicit expression for the function Sij, and the parametrization of p4 in terms of 
xi and X2 for a given collinear direction, it is easy to verify that when the collinear direction 
is the direction of the hard gluon gh, h G (1,2,3), terms in Eq. (4.26) that are proportional 
to \ contribute to Fi, terms proportional to X^-^^ 4 or contribute to F2 and terms that 

are proportional to F^j^^ with i ^ h, j ^ h, contribute to F3. 

As the next step, we check that the parametrization in Eq. (4.22) is consistent with the 
behavior of the real- virtual matrix elements in the collinear limit. Consider for definiteness 
the Sc^^^^ sector, where the singularity occurs when gluon g^ becomes collinear to gluon g^. 
In the collinear limit, color-ordered matrix elements factorize as follows [47] 



KL-^Hg..,. = 5^Split|- ,3.4 ® KL-^Hg. + ^^'Sphti;T939. 



M, 



(0) 



gig2-fHg3^ 



(4.27) 
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where the convolution sign refers to a sum over the hehcities of the intermediate gluon g^. 
The tree sphtting function for gp — )• gagt reads [47] 

Split^'''"' {gp gagb, z) = • ebh ■<^p + h- ^a^p -eb-ka- e^ea ■ ep) , (4.28) 

Sab 

where by definition z = Ea/{Ea + Eb) is the momentum fraction carried by ga, ^a,b,p are the 
polarization vector of three particles that participate in the splitting and all momenta are 
taken to be outgoing, so that p + ka + kb = 0- The one-loop splitting function is given by 



Spliti-'°°P(5p ^ gagb, z) = ^ Ftrec(^) Split*-^(gp ^ gagb, z) 
+ r- ^ F^cw{z) {ka - kb) ■ ei ( Sab^a ' ^b '^kb " (-aka ' ^b) ■ 

In conventional dimensional regularization, the functions -Ftrec and -Fnew read 

Ftree = [zfl{z) + (1 " z)h{l - z) - 2h] , 

2 \-SabJ 

^-w- (l_2e)(3-26) [sabj 



(4.29) 



(4.30) 



where 



fi = \cr f-r(l - e)r(l + e)z-i-^(l - zf -- + ^p^^e, e, 1 + e, z] 

\ z z 

/2 = --cr. 



(4.31) 



In calculating the splitting functions, one has to be careful with imaginary parts. Note that 
Sab ii^ Eq.(4.30) can be both positive and negative so that {—Sab)~'' may or may not have 
an imaginary part. Also, because of analytic continuation in the space-like region z can be 
smaller or larger than one, in which case /i may have an imaginary part. However, these 
cases are mutually exclusive since Sab is positive for the final state splitting where z < 1, and 
negative for the initial state splitting where z > 1. As a result, we do not need to care about 
the interference of two imaginary parts, or about their consistent definition. We require the 
expansion of the hypergeometric function in Eq.(4.31) through 0{e'^). It reads 



2Fi{e,e,l + e,z) = 1 + U2{z)e^ + e' 

-Lh{l -z)-U3{z) 



Cs + ^lnzln2(l - z) +ln(l - z)Li2(l - z) 



(4.32) 

+ 0{e^). 



We will need products of splitting amplitudes summed over polarization states of unre- 
solved particles. These polarization states must be taken in d-dimensions but, because of the 
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real-virtual kinematics, the four-momenta of all gluons are four-dimensional. We write these 
products as 



Splii'^'^'ig^ ^ 5a5fe)Split*'-°'=(gp^ ^ gaOb) = —P^^^^'^^'^iz, k,, e), (4.33) 



AajAj, ''''' 



with the LO splitting function defined in Eq. (4.6) and 

Split'^^'ig^ ^ 5a56)Spliti-i°°P(<7p^ ^ ga9b) = —P^g\^,{z,Ka,e), 
Xa,Xt """^ (4.34) 

With these definitions, we are in position to present the limiting behavior of the inter- 
ference of one-loop and tree amplitudes in the collinear limit. We find 



"^^^ y^^gg-^Hgg-^^^gg^Hgg) " W"^^ t,,gg^Hg-^^^u,gg^Hg ) 



Note that in the last term we have taken the interference splitting function outside of the 
real part and have replaced this splitting function by its real part. We are allowed to do that 
because Pggi^t is a symmetric tensor, so that we can write 



■^^^f,,gg^Hg-^^^u,gg^Hg^gg,mt " 2 f^,gg^Hg-^^^u,gg^Hg + l^* ^ ^))^gg,int 



y^^fi,gg^Hg-^^^u,gg^Hg) ^gg,inf 



(4.36) 



This observation is useful since we need Pgg i^x. both smaller and larger than one, and the 
above equation implies that the analytic continuation of Pgg\T^yt{z) can be done in an arbitrary 
way since the imaginary part drops out. 

It follows from Eq. (4.35) that in the collinear limit the amplitude has a single branch cut 
Indeed, Pgg{z) is a rational polynomial of z and therefore contributes to Fi while Pgg i^t 
is proportional to ~ x^*^. We therefore match the X2 — >• limit to Fi and F2 and require 
that ^3(2:1,0) vanishes. Finally, we note that the splitting functions in Eq. (4.35) exhibit 
spin correlations; we can handle them in exactly the same way as described in the Section 
dedicated to next-to- leading order computations, where we explained that for each phase- 



space point we compute the vector such that kj_ = y —k'^K^. We do this at the level of 
phase-space point generation, where we resolve all singularities related to the collinear fc^ — t- 
limit analytically. Once the vectors are known, we can rewrite Eq. (4.35) through sums 
over (d-dimensional) helicities in complete analogy with what was done at next-to-leading 
order. 
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5 Higher-order e terms in amplitudes 

The computational algorithm that we discuss in this paper is based on conventional di- 
mensional regularization, in which the polarization states of all particles are continued to 
d-dimensions. Therefore, it becomes an important issue in our construction to understand 
how scattering-amplitude contributions at higher orders in e can be calculated. The goal of 
this Section is to discuss this issue. 

We begin by pointing out that the highest-multiplicity amplitudes at any order of per- 
turbation theory are needed to 0{e^), since they only contribute to the finite parts of the 
relevant correction. Therefore, at NNLO for Higgs plus jet production, we can use four- 
dimensional expressions for tree amplitudes gg — )• Hggg and we can truncate one-loop am- 
plitudes gg — )• Hgg at 0{e^). However, for lower-multiplicity amplitudes, such as tree-level 
gg — )• Hg and gg — t- Hgg, we need to know higher-order e terms. In principle, it can be 
expected that higher-order e terms for the one- loop gg — )• Hg amplitudes are needed but, as 
it was pointed out in Ref. [68], this is not the case. Indeed, the 0{e) contributions to these 
amplitudes cancel out between the two-loop virtual correction, the square of the one-loop 
amplitude and the singular limit of the real-virtual correction. We use this cancellation as 
a consistency check on our numerical implementation. We calculate the one-loop gg — t- Hg 
amplitude through 0{e^) and check that the higher-order e terms do not contribute to the 
final result due to the above-mentioned cancellations. 

We note that computations of matrix elements squared for d 7^ 4 are straightforward 
if they are performed by adding and squaring Feynman diagrams. All one needs to do in 
this case is to use the correct contractions of metric tensors obtained after summing over 
polarization states of external particles. Unfortunately, if a calculation is done in this way, 
the results rapidly become unwieldy, especially when a large number of gluons is involved. 
Instead, we decided to compute higher-order e terms directly at the amplitude level. This is 
possible because for lower-multiplicity final states we can choose to parametrize momenta of 
all particles as four-dimensional. The "extra-dimensional" polarization vectors then have a 
simple property that they are orthogonal to all momenta, eg-pi = 0. Therefore, the only way 
such polarizations can contribute to the amplitude and give non-vanishing contributions is 
through scalar products ei^s ■ ej^s' = ~^s,s'- This implies that the necessary condition for the 
amplitude to be non-vanishing is that an even number of particles has "extra-dimensional" 
polarizations. 

To illustrate how this works in detail, we consider the tree-level gg — t- Hg amplitudes. 
We write the full amplitude as 

^(l'^^2^^3'^^) = 2ag!,^5s(F^^)c,c3^(l'^2"^3'^^), (5.1) 

where A^^^ is the Higgs effective coupling defined in Eq. (2.4) and (F'^2)ciC3 = —iV^f'^^'^^'^'-^ is 
the color generator in the adjoint representation. Apart from normal helicities, we can have 
amplitudes where exactly one pair of gluons has identical extra-dimensional polarizations. 
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We find one such independent color-ordered amplitude which reads 

A(l^2^3+ 



[13] [32] 2 (12) 



[21] 



(23) (31) 



(5.2) 



Because extra-dimensional polarizations in the above amplitude should be the same and 
because there are d — A = — 2e extra-dimensional directions, the matrix element squared for 
gg — >• Hg can be written as 

2 



2e |^(l'',2^3*)|V |^(1^2'^,3^)|V |^(1^2^3'^)|^j 



(5.3) 



=(l-e) \A{l^\2^\'i^--') 



The final result here can be easily verified since it implies that, except for the — 4m| term, 
the 0{e) contribution to the squared matrix element for gg — )• Hg and the 0{e^) contribution 
coincide up to a sign. 

We note that, in addition to the matrix element squared, our construction requires more 
complicated objects that appear in coUinear limits 



iA^spin(n)|2 = ^(l^2'^^3'^3) A*(l^2^^3'^3), 

h2,h3 



(5.4) 



where the amplitudes on the right-hand side are computed under the assumption that the 
polarization vector of the gluon gi is n^. To calculate |A4spin(f^)|^, we write n as a linear 
combination of suitable polarization vectors of gluon gi. In doing so, it is important to 
remember that the vector n can have extra-dimensional components. We write 



|A^spm(n)|^ = Y pin,hi,h[)\Mspin{hi,h[)\'^, 

hi,h[ 

\Mspin(.h,h[)\^ = J2 ^(l'^^2'^^3'^^) ^*(l'^^2'^^3'^^), 



(5.5) 



where p{n,h,h') = (n • Ch) {n ■ C/j/). The helicity labels hi,h'^ can assume the following 
values: {hi,h[) = (ij), (i, s), {s, j), {s, s) where i,j = ±. It is straightforward to compute 
\-Mspm{hi, h[)\'^ for all pairs of helicity labels. The key point here is that, for non-vanishing 
amplitudes, there must be either zero or two s-helicity labels. We are then left with the 
following non-zero entries 



|Mspin(i, j) 



|-^spin|£;_4 ih 



j)-2e [A{l\2',3')A*{l^,2',3' 



|M,pin(s,5)l' = 5^[^(l^2^3'^)^*(l^2^3^)+^(l^2^3^)x(l^2^3^)] . 



(5.6) 
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Finally, we note that following a similar approach, it is straightforward to obtain the double- 
correlated matrix element |Mspin(/ii, /i'^ , /i2) ^2)^' which is needed to describe singular limits 
in the double-collinear sectors. 

We also need to discuss the e-dependent parts of — > Hgggg amplitudes. In this case, 
we use the following color decomposition 

(TG52 

The situation now is slightly more involved than before because there are more options for 
extra-dimensional polarizations. Indeed, with four gluons the amplitude does not vanish 
if all of them have identical extra-dimensional polarizations but also when there are two 
pairs of gluons with different extra-dimensional polarizations. We will denote color-ordered 
amplitudes for these cases as 2*, 3*, 4*) and A(l'*,2*,3'*',4*'). These amplitudes can be 
written in a relatively compact form. For example, 

4 

2^ 3^ 4*) = ^ n'F{l, 2, 3, 4), 

i=0 

F(l,2,3,4) = ^ 1 + ^ + ^ - -^ + -^ + , 

S123 V ■523 ■512/ \^Si2 2S23j ^ V ^14523 S12S34/ 

2^ 3^', 4^') = ml ( £ii±ii^ - -ii^ + ^ii±^ - --^) ^^-^^ 



S12S124 ■Sl2'Sl23 ■S34^Sl34 ■S34S234, 
3£24 

.Sl2^S34 .S12S34 



_^ gl4.S23 _ ■S13S24 _ ^ 



A{r,2'',3',A'') = 2-ml [— + — + — + 



SUi S134 S234 ■5123 

where 7^ is a permutation operator defined as TlF{a,b,c,d) = F{b,c,d,a). The amplitudes 
remain compact even if only one pair of gluons has extra-dimensional polarization. For 
example, we obtain 

. s OS 4+^ ^ {l\ph\mPh\M , {l\ph\2]mph\2] mlilSf 
^ ' ' si23(12)(23) + ^i34(14)(34) + (12)(14)(23)(34)' ^^'^^ 

where ph is the outgoing momentum of the Higgs boson. Similar results for all other helicity 
configurations can be derived. 

We are now in position to discuss how to use these amplitudes to assemble the matrix 
element squared for — ?• Hgggg, summed over polarization vectors of all gluons. Similar 
to the —7- Hggg case that we already discussed, amplitudes with two gluons with extra- 
dimensional polarizations, e.g. ^(1*, 2-^, 3^,4'*), enter with a (d — 4) = — 2e weight. The same 
is true for the amplitude 2*, 3*, 4*), as s just counts the number of extra-dimensional 

polarizations. For amplitudes like ^(1*, 2*, 3*', 4*'), we have again d — 4 polarizations for the 
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index s and d — 5 for s' since, by construction, s ^ s' . Combining everything, we obtain 



\M{H,gi,g2,g3,gA)\'^ = \M{H,gi,g2,g3,gA)\l=A - 2e 



l^(l^2^3^4 



sm2 



+ (|^(l^2^3'^^4'^oP + l-^(l^2''^3^4'^OI^ + l-^(l^2''^3''^4 



sm2 



(5.10) 



hi as\\2 



+ 2e(2e + l) 2^ 3"', 4''')|2 + |^(1^ 2"' , 3^ 4"' )|^ + 2*' , 3"' , 4 



s'\\2 



+ 



is' AS\\2 



In full analogy with the amplitudes for — )• Hggg, we can calculate \M{hi, /i'i)spinP which 
is required to describe spin correlations in collinear limits. A simple analysis reveals that this 
spin-correlated matrix element squared is non- vanishing provided that [hi,h'^) = or 
{hi,h'^) = (s, s) so that no mixed terms as (i, s) or {s,s') appear. The result can be written 
as 

iM(z,j%pinp = iM(z,j%piniL4 - 26 5^ (^(l^2^3^4V*(l^2^3^4^) 

h=± 

+A{i\ 2^ 3^ 4^)X(i^ 2^ 3^ 4^) + a{i\ 2\ 3^ 4^)X(i^ 2\ 3^ 4^) 



|M(s,s) 



spml 



2^ 3^ 4'^)x(l^ 2^ 3^ 4^) + ^ 2*, 3^ 4^)^*(l^ 2^ 3^ 4^ 



+A{V, 2\ 3^ A^)A*{V, 2\ 3^ 4^) + A{V, 2^ 3\ A^)A*{1' ,2' ,3\A^)) 

- (i + 2e) (^(l^2^3*',4^')^*(l^2^3''',4^')+^(l^2*',3^4*')x(l^2''',3^4''') 

2'' , 3^' , 4^)^* (l^ 2^' , 3'*' , 4^)) . 
6 Numerical implementation 



(5.11) 



In this Section we discuss the implementation of the algorithm described above in a numerical 
program. We choose to do so in FORTRAN 90 since it offers the option of performing 
computations in double- and quadruple precision in a straightforward way. Such flexibility is 
important because in our framework singular limits are approached numerically, and we have 
to find a balance between the speed of the code and the numerical stability which requires 
switching to quadruple precision computations when close to singularities. 

For numerical implementation of the required amplitudes we used, as much as possi- 
ble, pieces of the FORTRAN 77 code MCFM [66]. After translating to FORTRAN 90 we 
checked our numerical implementation of the tree-level amplitudes for — t- Hggg, — t- Hgggg 
and —7- Hggggg processes against MadGraph [69]. As we explained earlier, since we work 
in conventional dimensional regularization, we need to know 0(e) parts of tree-level ampli- 
tudes, which are presented in the previous Section. These 0(e) parts were checked against a 
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Feynman diagram-based computation of amplitudes squared where explicit sums over gluon 
polarizations were performed. The relevant diagrams for — )• + ng, n = 3, 4 were obtained 
with QGRAF [70] and manipulated with FORM [71]. Finally, we note that since we require 
the one- loop corrections to gg — t- Hg through 0{e'^), we recomputed the one- loop gg — t- Hg 
amplitudes and compared them against the results presented in [72]. For the — )■ Hggg 
one-loop amplitudes, we borrowed significant parts of the FORTRAN code from MCFM. The 
one-loop integrals that are required for this calculation are computed using QCDloops [73]. 
The box one-loop master integral gg — t- Hg is needed to higher orders in the expansion in e, 
and can be obtained starting from an all-orders result in Ref. [38]. 

A central part of the described computational algorithm is the calculation of integrals of 
the following form 

1 

J dxi. . dxn dyT>i-^{xi)... T>i^{xn)F{xi, ..Xn,y), (6.1) 


where n counts the number of singular phase-space variables ( n = 4 for double-real and 
n = 2 for real- virtual) , y collectively denotes all non-singular variables, the functions P(x) 
are defined as . ^ 

Vo{x) = 6{x), V,{x) = , (6.2) 

and "^jij ^ 3. The function F(xi, ..Xn,y) is obtained by multiplying the matrix element 
squared for a particular physics process by appropriate powers of xi,..a;„, as explained in 
Section 3. To compute multi-dimensional integrals of the type shown in Eq. (6.1), we use the 
adaptive Monte-Carlo algorithm VEGAS [74] as implemented in the CUBA library [75]. We 
note that when plus distributions are expanded out in Eq. (6.1), we obtain integrands that 
are iterations of the following basic form 



X 



^ [F{xi, . . . Xi-i,Xi, Xi+i ...) - F{xi, . . . 0, Xi+i, ...)]. (6.3) 



To understand subtleties of the numerical implementation of Eq. (6.1), it is important to 
realize that the two terms in the numerator of Eq. (6.3) are computed differently in the 
numerical code. Indeed, the function F{xi, . . . Xj, Xj+i . . . ) is calculated from the ma- 
trix element squared that describes the highest-multiplicity process for a given channel (for 
example, it is — )■ Hggggg for the double-real emission processes). On the other hand, 
F{xi, . . . Xj-i, 0, Xj+i . . . ) is computed by first analytically calculating the appropriate singu- 
lar limit from the full matrix element and then implementing that limit as an independent 
function or subroutine in the numerical code. This implies two things. First, 

hm F(..Xi,...) =F(..,0,..), (6.4) 

is an important and non-trivial check of the calculation and of its implementation in the 
numerical program. Second, because the full matrix elements become numerically unstable 
for very small values of x, it is not possible to calculate integrands as in Eq. (6.1) all the 
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way to xi^„^n = 0. In our numerical implementation, we follow the approach of Ref. [56] and 
require that the product of all generated singular variables is larger than a small parameter 



Independence of the final result from the value of 6c is an essential check of the correctness 
of the numerical implementation; we will discuss the necessary condition for that in the next 
Section. 

To obtain results for partonic cross-sections that will be presented in the next Section, 
we use 5c = 10"^'^. We introduce a switch in the program that forces a quadruple precision 
calculation of the integrand in Eq. (6.1) to occur provided that xiX2--Xn < 6s, where 5s is 
conservatively chosen to be 5s = 10"^. We find that, compared to a pure double-precision 
computation, our implementation of the switch slows a calculation by about a factor of two. 
This, however, is not a problem since the program is quite fast because it employs helicity 
amplitudes to construct the relevant matrix elements. To illustrate how fast the program is, 
we note that to get contributions to integrated partonic cross-sections for one center-of-mass 
collision energy, we need about half an hour to obtain all the poles in e and about four hours 
to obtain all the relevant finite parts provided that a calculation is done on a cluster of twenty 
eight-core 2.83 GHz nodes. 

Finally, we note that our current implementation of the numerical integration procedure 
allows us to calculate partonic cross-sections but not kinematic distributions. This, however, 
seems a relatively minor problem since at every step of the calculation we know the kinematics 
of the final state and we can access the weight. It appears therefore that the current imple- 
mentation can be easily extended to make a true parton-level generator capable of computing 
different observables in a single run. In fact, the possibility to do this within the current 
framework was recently demonstrated for the simpler processes t — t- be'^u and b — ?• uev' in 
Refs. [76, 77]. We plan to return to the discussion of this issue in the context of Higgs boson 
production in the near future. 

7 Checks and final results 

In this Section, we describe checks on the calculation and present the results for the partonic 
cross-section. The first check that we describe follows from the fact that, in the numerical 
program, Eq. (6.4) is non-trivial to satisfy, because of the different ways in which the function 
F{xi,X2, ..) and its boundary values are computed. On the other hand, Eq. (6.4) is a necessary 
requirement for the existence of integrals shown in Eq. (6.1), so that its validity in our 
numerical program should be carefully investigated. To check Eq. (6.4), we compute 



as a function of t — )• 0, for random choices of x = [xi,X2,---]- If the calculations are done 
properly, we should find Lij {t) — )• as t — )• 0, independent of x. Also, because the variables 



XiX2..Xn > 5c. 



(6.5) 




F(^Xl, . .tXi^ , .. , tXi2 , . .) 
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X define the kinematics of the process, each function Li^^i^, {t) probes a particular singular 
limit of the full amplitude. 

For the double real emission sectors, we consider fifteen different limits, for example Lx\ 5 
Lx2, Lx-j--, Lxi^x2---j Lxi^x3,x4,, and check numerically how these functions approach zero. In 
particular, we know that all the soft limits should scale as t, while collinear limits should scale 
as \/t. To illustrate this point we plot distributions for the functions Lx^it) and Lx^{t) in 
Fig. 1, for two sample sectors. The function Lx^{t) describes the soft limit and the function 
Ljjg describes the collinear limit. To obtain these plots, ten thousand x points were randomly 
generated and the two functions Lx^{t) and Lx^{t) were computed for two values of t that 
differ either by one (soft) or two (collinear) orders of magnitude. It is evident from Fig. 1 that 
the widths of the resulting distributions scales with the parameter t as expected. We also 
note that, in case of the collinear limit, the quality of the distribution is very sensitive to the 
correct implementation of spin correlations. In fact, by removing the spin-correlation part 
from collinear splitting functions, we find Lx^{t) ~ 0(10"^) independent of t for t < 10~^. 



Double soft limit, 4||1 , 5\\2 Triple collinear limit, 4||3, 5||3, sector 1 




L(x,) L(X3) 



Figure 1. Scaling behavior for soft (left) and collinear (right) double-real emission limits, as obtained 
with our Fortran code in quadrupole precision. See the text for explanation. 

We note that we cannot follow the same strategy to check the 0{e) terms for lower- 
multiplicity amplitudes, since we do not have a computation of the — )• Hggggg amplitude 
beyond 0{e^). We can nevertheless check the consistency of our calculation and implemen- 
tation by comparing different limits against each other. In total, we consider 60 different 
combinations for all double-real sectors and check that each of them behaves in a way that is 
similar to what is shown in Fig. 1, for e = 0, 1, 2. 

To check the implementation of the real- virtual corrections, we need to modify the above 
strategy, since Fiiy{x) is given by a linear combination of three functions with potentially 
logarithmically-singular coefficients, as shown in Eq. (4.22). To probe soft and collinear limits 
in the real-virtual case, we define two functions 

r / .X 1 T; [FRv(tei,X2,...)] r / .^ 1 Te [FRv{e,xx,tx2,...)] , . 
Li{e,t} = 1 - -, L2[e,t) = 1 - — , (7.2) 

7e Gi(t,xi,X2, ...) Te G2 (e, t, xi , a;2 , ..) 
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where 

Gi{t,Xi,X2, ..) 
G2{t,Xi,X2, ..) 



Fi{0,X2,...)+F2{0,X2,. 

Fi{xi,0,...)+F2{xi,0,. 



+ F3(0,X2,...) [t^x^^ 



(7.3) 



XiX2\ 
[txlx2] 

The operator Te in Eq. (7.2) imphes that the relevant term in the Laurent expansion in e of 
the corresponding function should be taken. For illustrative purposes, we show distributions 
of Lxj^{0,t) and Lx2iO,t) for one of the sectors in Fig. 2. Similar to the double real emission 
case, we observe the 0{t) scaling in the soft limit and the 0{^/t) scaling of the collinear limit. 



Real virtual sector 4| [3 



Real virtual sector 4||1 



700 





1x10"' 



Figure 2. Scaling behavior for soft (left) and collinear (right) real- virtual emission limits, as obtained 
with our Fortran code in quadrupolc precision. See the text for explanation. 

A further check of the correctness of the calculation is provided by the cancellation 
of poles. Singularities of double-real, real-virtual and double-virtual contributions start at 
0{e~'^). Starting from order 0(e~^), collinear subtractions, renormalization and contributions 
related to extra-dimensional components of the unresolved momenta are required for the 
cancellation of poles. We note that within our framework, we compute coefficients of the 
Laurent expansion in e and check the cancellation of poles numerically. To see how well this 
cancellation works, we compute the ratios 

c o-RR(e) + o"Rv(e) + crvv{e) + (e) -I- CTrenorm(e) + 0"d-4(e) .x 

|o-RR(e)| + |(TRv(e)| + |<Tyy(e)| + | 

•^conv 

(e)| + |(Trenorm(e)| + 1 0-^-4 (e)| 

at various orders in e. In Eq. (7.4), we account for double-real, double- virtual, real- virtual 
contributions as well as convolutions, renormalization and the contribution due to extra- 
dimensional components of the unresolved gluon momenta. We show 5{e) in Fig. 3 for e = —2 
and e = — 1. Interestingly, it appears from Fig. 3 that we loose almost one order of magnitude 
in the quality of cancellation when we move from 0{€~^) to 0{e~^). Nevertheless, at 0{e-^) 
the cancellation is at the level of few per mille or better, which is acceptable. Finally, we 
note that omission of extra-dimensional components in the momentum parametrization leads 
to residual non-cancellation of singularities at the level of (52 ~ 5 x 10~^ and 5i ~ 2 x 10~^, 
which is very large compared to values of 6 that we observe in Fig. 3. 
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Figure 3. Residuals of poles in e for the total cross-section as the function of partonic center-of-mass 
energy. The left panel shows ©(e^^), and the right panel shows 0{e~^). See the text for explanation. 



As a final check of the calculation, we discuss the dependence of the result on the renor- 
malization and factorization scales. In this paper, we equate them and denote both by /x. 
We can compute the /i-dependence of the cross-section either by introducing /u^ per cou- 
pling constant in the various elements of the calculation in the standard way, or by solving 
the renormalization group equation that follows from the fact that convolution of the par- 
tonic cross-section with parton distribution functions is /i-independent. The results of this 
computation can be found in Section 2. We have checked that when the /i-dependence is 
computed with our numerical code, the result agrees with the analytic computation based on 
renormalization group invariance. 

We now present our results. We compute the hadronic cross-section for the production of 
the Higgs boson in association with a jet at the 8 TeV LHC through NNLO in perturbative 
QCD. We reconstruct jets using the /c_|_-algorithm with AR = 0.5 and p±j = 30 GeV. The 
Higgs mass is taken to be niH = 125 GeV and the top-quark mass nit = 172 GeV. We use 
the latest NNPDF parton distributions [78, 79] and numerical values of the strong coupling 
constant at various orders in QCD perturbation theory as provided by the NNPDF fit. We 
note that in this case as{mz) = [0.130,0.118,0.118] at leading, next-to-leading and next-to- 
next-to-leading order, respectively. We choose the central renormalization and factorization 
scales to he = = n^H- In Fig- 4 we show the partonic cross section for gg ^ H + j 
multiplied by the gluon luminosity through NNLO in perturbative QCD 



^dohad _ ^ da{s,as,iJ.R,iJ.F) ^ ^ 



Shad 



(7.5) 



where /3 measures the distance from the partonic threshold. 



/3 = y 1 - Eth = ^Jml +pI^. +p^^^ ^ 158.55 GeV. (7.6) 

The partonic luminosity C is given by the integral of the product of two gluon distribution 
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Figure 4. Results for the product of partonic cross-sections gg ^ H + jet and parton luminosity in 
consecutive orders in perturbative QCD at fiR = fip = ruh = 125 GeV. See the text for explanation. 

functions 

dx f z \ 

C{z,^iF)= — 9ix,fJ-F)g [-,fJ'F) ■ (7.7) 

It follows from Fig. 4 that NNLO QCD corrections are significant in the region ^/s < 
500 GeV. In particular, close to partonic threshold ^/s ~ Eth, radiative corrections are en- 
hanced by threshold logarithms In /? that originate from the incomplete cancellation of virtual 
and real corrections. There seems to be no significant enhancement of these corrections at 
higher energies, where the NNLO QCD prediction for the partonic cross-section becomes al- 
most indistinguishable from the NLO QCD one. Note that we extend the calculation of the 
NNLO partonic cross-section to ^/s ~ 500 GeV only. From leading and next-to-leading order 
computations, we know that by omitting the region ^/s > 500 GeV, we underestimate the 
total cross-section by about 3%. To account for this in the NNLO hadronic cross-section cal- 
culation, we perform an extrapolation to higher energies constructed in such a way that when 
the same procedure is applied to LO and NLO cross-sections, it gives results that agree well 
with the calculation without extrapolation. The correction for the extrapolation is included 
in the NNLO QCD cross-sections results shown below. 

We now show the integrated hadronic cross-sections for the production of the Higgs 
boson in association with a jet at 8 TeV LHC in the all-gluon channel. We choose to vary 
the renormalization and factorization scale in the range = fip = ttih/'^, trtH, 2m h- After 
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Figure 5. Scale dependence of the hadronic cross section in consecutive orders in perturbative QCD. 
See the text for details. 

convolution with the parton luminositites, we obtain^ 

(TLo(pp^^i) = 2713t^2i6 fb, 

<TNLo(pp ^ Hj) = 4377;^60 ^7_g) 

(^NNLoiPP Hj) = 6177:^242 ft). 

We note that NNLO corrections are sizable, as expected from the large NLO -fC— factor, but 
the perturbative expansion shows marginal convergence. We also evaluated PDFs error using 
the full set of NNPDF replicas, and found it to be of order 5% at LO, and of order 1-2% at 
both NLO and NNLO, similarly to the inclusive Higgs case [78]. The cross-section increases 
by about sixty percent when we move from LO to NLO and by thirty percent when we move 
from NLO to NNLO. It is also clear that by accounting for the NNLO QCD corrections we 
reduce the dependence on the renormalization and factorization scales in a significant way. 
The scale variation of the result decreases from almost 50% at LO, to 20% at NLO, to less 
than 5% at NNLO. We also note that a perturbatively-stable result is obtained for the scale 
choice /u ~ m/f/2. In this case the ratio of the NNLO over the LO cross-section is just 1.5, 
to be compared with 2.3 for fi = mn and 3.06 for /u = 2mH, and the ratio of NNLO to NLO 
is 1.2. It is interesting to point out that a similar trend was observed in the calculation of 
higher-order QCD corrections to the Higgs boson production cross-section in gluon fusion. It 
has been pointed out that because of the rapid fall of the gluon PDFs, the production cross 
section is dominated by the threshold region, thus making /i = m/f/2 an excellent choice for 
the renormalization and factorization scales [14, 81]. The reduced scale dependence is also 
apparent from Fig. 5, where we plot total cross-section as a function of the renormalization 
and factorization scale /i in the region p_ij < n < 2mh- 

^We checked our LO and NLO resuhs against MCFM (gluons only), and found agreement. 
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Finally, we comment on the phenomenological relevance of the "gluons-only" results for 
cross-sections and X-factors that we reported in this paper. We note that at leading and next- 
to-leading order, quark-gluon collisions increase the H + j production cross-section by about 
30 percent, for the input parameters that we use in this paper. At the same time, the NLO 
X-factors for the full H+j cross-section are smaller hy about 10—15 percent than the "gluons- 
only" X-factors, presumably because quark color charges are smaller than the gluon ones. 
Therefore, we conclude that the gluon-only results can be used for reliable phenomenological 
estimates of perturbative X-factors but adding quark channels will be essential for achieving 
precise results for the H + j cross-section. We plan to return to this issue in the future. 

8 Conclusions 

In this paper we reported a calculation of the NNLO QCD corrections to the partonic process 
gg —7- H+jet. This is one of the first calculations where NNLO QCD corrections are computed 
to a 2 — 7- 2 process whose cross-section depends on the implementation of the jet algorithm 
already at leading order. We believe that gg — )• Hg is a sufficiently typical process to expose 
all non-trivial features of a generic NNLO computation for a 2 — ?■ 2 process at a hadron 
collider. Indeed, we have used this process to show that the computational technique that we 
describe in this paper can successfully deal with: 

• a large number of contributing Feynman diagrams; 

• colored particles in the initial and in the final state; 

• collinear subtractions and parton distribution functions; 

• all soft and collinear limits; 

• known helicity matrix elements; 

• spin correlations; 

• a realistic jet algorithm. 

The only "non-generic" feature that we benefited from is a much simpler bookkeeping that 
is required for gg Hg compared to the general case computation. 

We believe that the techniques reported in this paper that built upon earlier work de- 
scribed in Refs.[55, 56, 63], allow computation of the NNLO QCD corrections to an arbitrary 
2—7-2 process at hadron colliders provided that the corresponding two-loop matrix elements 
are available. Since this is the case for most of the processes that are desirable to know at 
NNLO (c.f. the "NNLO wishlist" in Ref. [80]), our results open up a way to perform the 
required calculations. 

On the other hand, it is not entirely clear to us how to extend the computational tech- 
nology reported in this paper to make it practically applicable to 2 — >• n, n > 2 processes. 
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In this case, the problem is related to the 0(e) parts of the amplitudes and the choice of 
extra-dimensional components to parametrize four-momenta of unresolved gluons. The point 
is that in the 2—7-2 process these details can still be dealt with by brute force, as we did 
in this paper, but for large n this will be increasingly difficult to do. Therefore, it is an 
interesting theoretical question to re-formulate this technique in such a way that much of the 
irrelevant 0(e) dependencies is avoided. We hope to return to this point in the future. 
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A Appendix 

We report here the formulae for the splitting functions and their convolution needed for the 
renormalization of parton distribution functions at NNLO, as described in Sec. 2. 



11 
12 



S{l-x) + 



22 

y 



1 



1 — X 



1 — X 



+ 



N 1 — X 

+ x(l - X H 1 

X 

ln(l /121 2 2,,,, 

— ^ +{ vr^ U(l-x) 

l-x , \ 36 3 ' ^ ' 



+ 



4{x^ + 3x^ + 1- 4x^ 



x{x — 1) 



■ ln(x) 



8(2x - 1 + - 



67 



■ln(l - x) + 



2(llx3-7x2-4x-ll) 



3x 



- + 3C3)5(l-x)+( -^-2C2 



1 — X 



4(x2+x + l)2 

+ TT——\ Li2(-x) 

x(l -|- X) 



+ 



+ 



4(3 2x2 + 4x 2x3) ^ ^4(x2-x-l)2 



2(1 + x) 
4(x2 J^x + Vf 



<2 + 



2(1 -x2) 



- + 


x( 


25 


109 
1^" 


18 ~ 



, , , 4(x2 - a; + 1)2 75 33 44x2 . ^ ^ 

mil + x) ; ^ ln(l — x) 1 X In(x) 

x(l + x) ^ ' x(l-x) ^ ^ 9 9 3 ' ^ ^ 



(A.l) 



References 

[1] G. Aad et al. [ATLAS Collaboration], Phys. Lett. B 716, 1 (2012) [arXiv:1207.7214 [hep-ex]]. 



-43- 



[2] S. Chatrchyan et al. [CMS Collaboration], Phys. Lett. B 716, 30 (2012) [arXiv:1207.7235 
[hep-ex]]. 

S. Chatrchyan et al. [CMS Collaboration], arXiv:1212.6639 [hep-ex]. 

LHC Higgs Cross Section Working Group, A. David, A. Denner, M. Duehrssen, M. Grazzini, 

C. Grojean, G. Passarino and M. Schumacher et ai, arXiv:1209.0040 [hep-ph]. 

S. Dawson, Nucl. Phys. B359, 283 (1991). 

A. Djouadi, M. Spira and P. Zerwas, Phys. Lett. B264, 440 (1991). 

D. de Florian, M. Grazzini and Z. Kunszt Phys. Rev. Lett. 82, 5209 (1999). 
V. Ravindran, J. Smith and V.L. van Neerven, Nucl. Phys. B634, 247 (2002). 
C. J. Glosser and C.R. Schmidt, JHEP 0212, 016 (2002). 
J.M. Campbell, R.K. Ellis and G. Zanderighi, JHEP 0610, 028 (2006). 
J. M. Campbell, R. K. Ellis and C. Williams, Phys. Rev. D 81, 074023 (2010). 
H. van Deurzen, N. Greiner et al., arXiv: 1301. 0493 [hep-ph]. 

R. V. Harlander and W. B. Kilgore, Phys. Rev. Lett. 88, 201801 (2002) [hep-ph/0201206]. 

C. Anastasiou and K. Melnikov, Nucl. Phys. B 646, 220 (2002) [hep-ph/0207004]. 

V. Ravindran, J. Smith and W. L. van Neerven, Nucl. Phys. B 665, 325 (2003) 
[hep-ph/0302135]. 

C. Anastasiou, K. Melnikov and F. Petriello, Nucl. Phys. B724, 197 (2005). 

S. Catani and M. Grazzini, Phys. Rev. Lett. 98, 222002 (2007); M. Grazzini, JHEP 0802, 043 
(2008). 

A. Banfi, G.P. Salam and G. Zanderighi, JHEP 1206, 159 (2012). 
T. Becher and M. Neubert, JHEP 1207, 108 (2012) [arXiv;1205.3806 [hep-ph]]. 
F. J. Tackmann, J. R. Walsh and S. Zuberi, arXiv:1206.4312 [hep-ph]. 
X. Liu and F. Petriello, Phys. Rev. D 87, 014018 (2013) [arXiv:1210.1906 [hep-ph]]. 
T. Kinoshita, J. Math. Phys. 3, 650 (1962). 
T. D. Lee and M. Nauenberg, Phys. Rev. 133, 1549 (1964). 

S. Catani and M. H. Seymour, Nucl. Phys. B 485, 291 (1997) [Erratum-ibid. B 510, 503 
(1998)] [hep-ph/9605323]. 

S. Frixione, Z. Kunszt and A. Signer, Nucl. Phys. B 467, 399 (1996) [hep-ph/9512328]. 
K. Melnikov and F. Petriello, Phys. Rev. D74, 114017 (2006). 

K. Melnikov and F. Petriello, Phys. Rev. Lett. 96, 231803 (2006) [hep-ph/0603182]. 
S. Catani, L. Cieri, G. Ferrera, D. de Florian, M. Grazzini, Phys. Rev. Lett. 103, 082001 (2009). 
C. Anastasiou, K. Melnikov and F. Petriello, Phys. Rev. Lett. 93, 262002 (2004). 
C. Anastasiou, K. Melnikov and F. Petriello, Phys. Rev. Lett. 93, 032002 (2004). 



- 44 - 



[31] A. Gehrmann-De-Ridder, T. Gehrmann, E.W.N. Glover and G. Heinrich, JHEP 0712, 094 
(2007). 

[32] A. Gchrmann-Dc-Ridder, T. Gehrmann, E.W.N. Glover and G. Heinrich, JHEP 0711, 058 
(2007). 

[33] S. Weinzierl, Phys. Rev. Lett. 101, 162001 (2008). 

[34] C. Anastasiou, K. Mclnikov, F. Petriello, JHEP 0709, 014 (2007). 

[35] S. Biswas, K. Melnikov, JHEP 1002, 089 (2010). 

[36] K. Melnikov, Phys. Lett. B666, 336 (2008). 

[37] G. Ferrera, M. Grazzini and F. Traniontano, Phys. Rev. Lett. 107, 152003 (2011) 
[arXiv: 1107. 1164 [hep-ph]]. 

[38] C. Anastasiou, F. Herzog and A. Lazopoulos, JHEP 1203, 035 (2012) [arXiv:1110.2368 
[hep-ph]]. 

[39] S. Catani, L. Cieri, D. de Florian, G. Ferrera and M. Grazzini, Phys. Rev. Lett. 108, 072001 
(2012) [arXiv:1110.2375 [hep-ph]]. 

[40] S. Catani, D. de Florian, and M. Grazzini, Nucl. Phys. B596, 299 (2001). 

[41] F. A. Berends and W.T. Giele, Nucl. Phys. B313, 595 (1989). 

[42] J.M. Campbell and E.W.N. Glover, Nucl. Phys. B527, 264 (1998). 

[43] S. Catani and M. Grazzini, Phys. Lett. B446, 143 (1999). 

[44] Z. Bern, et al. Phys. Rev. D60, 116001 (1999). 

[45] S. Catani and M. Grazzini, Nucl. Phys. B570, 287 (2000). 

[46] S. Catani and M. Grazzini, Nucl. Phys. B591, 435 (2000). 

[47] D. A. Kosower, P. Uwer, Nucl. Phys. B563, 477 (1999). 

[48] C. Anastasiou, E.W.N. Glover, C. Olcari and M.E. Tejcda-Ycomans, Nucl. Phys. B 601, 318 
(2001) [hep-ph/0010212]; 601, 347 (2001) [hep-ph/0011094]; 605, 486 (2001) [hcp-ph/0101304]; 
E.W.N. Glover, C. Olcari and M.E. Tejcda-Ycomans, Nucl. Phys. 605, 467 (2001) 
[hep-ph/0102201]; 

C. Anastasiou, E.W.N. Glover and M.E. Tejcda-Ycomans, Nucl. Phys. B 629, 255 (2002) 
[hcp-ph/0201274]; 

E.W.N. Glover and M.E. Tejeda-Yeomans, JHEP 0306, 033 (2003) [hep-ph/0304169]; 
E.W.N. Glover, JHEP 0404, 021 (2004) [hep-ph/0401119]; 

Z. Bern, A. De Freitas and L.J. Dixon, JHEP 0109, 037 (2001) [hep-ph/0109078]; JHEP 0203, 
018 (2002) [hep-ph/0201161]; JHEP 0306, 028 (2003) [hep-ph/0304168]; 
A. De Freitas and Z. Bern, JHEP 0409, 039 (2004) [hcp-ph/0409007]. 

[49] L.W. Garland, T. Gehrmann, E.W.N. Glover, A. Koukoutsakis and E. Remiddi, Nucl. Phys. B 
627, 107 (2002) [hcp-ph/0112081] and Nucl. Phys. B642, 227 (2002) [hcp-ph/0206067]. 

[50] T. Gehrmann, M. Jaquier, E. W. N. Glover and A. Koukoutsakis, JHEP 1202, 056 (2012) 
[arXiv;1112.3554 [hep-ph]]. 

[51] K.G. Chetyrkin, B. Knichl and M. Steinhauser, Phys. Rev. Lett. 79, 353 (1997); Nucl. Phys. 
B510, 61 (1998). 



-45 - 



[52] S.D. Badger and E.W.N. Glover, P. Mastrolia and C. Williams, JHEP 1001, 036 (2010). 

[53] V. Del Duca, A. Frizzo and F. Maltoni, JHEP 0405, 064 (2004); L.J. Dixon, E.W.N. Glover, 
V.V. Khoze JHEP 0412, 015 (2004). 

[54] S. Catani and M. Grazzini, Nucl. Phys. B 570, 287 (2000). 

[55] M. Czakon, Pliys. Lett. B 693, 259 (2010). 

[56] M. Czakon, Nucl. Phys. B849, 250 (2011). 

[57] T. Binoth and G. Heinrich, Nuel. Phys. B585, 741 (2000). 

[58] T. Binoth and G. Heinrich, Nucl. Phys. B 693, 138 (2004). 

[59] C. Anastasiou, K. Melnikov and F. Petriello, Phys. Rev. D69, 076010 (2004). 

[60] P. Baernreuther, M. Czakon and A. Mitov, Phys. Rev. Lett. 109, 132001 (2012). 

[61] M. Czakon and A. Mitov, JHEP 1212, 054 (2012). 

[62] M. Czakon and A. Mitov, JHEP 1301, 080 (2013) [arXiv:1210.6832 [hep-ph]]. 

[63] R. Boughezal, K. Melnikov and F. Petriello, Phys. Rev. D 85, 034025 (2012) [arXiv:1111.7041 
[hep-ph]]. 

[64] A. Gehrmann-De Ridder, T. Gehrmann and E. W. N. Glover, JHEP 0509, 056 (2005) 
[hep-ph/0505111]; A. Daleo, T. Gehrmann and D. Maitre, JHEP 0704, 016 (2007) 
[hep-ph/0612257];J. Currie, E. W. N. Glover and S. Wells, arXiv:1301.4693 [hep-ph]; 
R. Boughezal, A. Gehrmann-De Ridder and M. Ritzmann, JHEP 1102, 098 (2011) 
[arXiv:1011.6631 [hep-ph]]; A. Gehrmann-De Ridder, T. Gehrmann and M. Ritzmann, JHEP 
1210, 047 (2012) [arXiv:1207.5779 [hep-ph]]. 

[65] A. G. -D. Ridder, T. Gehrmann, E. W. N. Glover and J. Pires, arXiv:1301.7310 [hep-ph]. 

[66] J. M. Campbell and R.K. Ellis, Phys. Rev. D 62, 114012 (2000). The MCFM program is 
publicly available from http:/ /mcfm. fnal.gov. 

[67] S. Catani and M. Grazzini, Nucl. Phys. B 591, 435 (2000). 

[68] S. Weinzierl, Phys. Rev. D 84, 074007 (2011). 

[69] J. Alwall, M. Herquet, F. Maltoni, O. Mattelaer and T. Stelzer, JHEP 1106, 128 (2011) 
[arXiv:1106.0522 [hep-ph]]. 

[70] P. Nogueira, J. Comput. Phys. 105, 279 (1993). 

[71] J. A. M. Vermaseren, math-ph/0010025. 

[72] C. R. Schmidt, Phys. Lett. B413, 391 (1997) [hep-ph/9707448]. 

[73] R. K. Ellis and G. Zanderighi, JHEP 0802, 002 (2008) [arXiv:0712.1851 [hep-ph]]. 

[74] G.P. Lepage, Cornell preprint CLNS-80/447. 

[75] T. Hahn, Comput.Phys.Commun. 168, 78 (2005). 

[76] M. Brucherseifer, F. Caola and K. Melnikov, arXiv:1301.7133 [hep-ph]. 

[77] M. Brucherseifer, F. Caola and K. Melnikov, arXiv:1302.0444 [hep-ph]. 



-46- 



[78] R. D. Ball, V. Bertone, S. Carrazza, C. S. Deans, L. Del Debbio, S. Forte, A. Guffanti and 
N. P. Hartland et ai, Nucl. Phys. B 867, 244 (2013) [arXiv:1207.1303 [hep-ph]]. 

[79] R. D. Ball et al. [NNPDF Collaboration], Nucl. Phys. B 855, 153 (2012) [arXiv: 1107.2652 
[hep-ph]]. 

[80] See a talk by J. Huston at the LoopFcst XI, Pittsburgh, May 2012. 

[81] M. Kramer, E. Laenen and M. Spira, Nucl. Phys. B 511, 523 (1998) [hep-ph/9611272]. 



-47- 



